= tee 
: =) ; 
— @ =z = 
a: w) = a 
ae = 2 
Sea @ 
@® Ee Ser 
Z 
as is 
= a 
-, 
Recy 2 
ea 
© A 
Ratbsss © 
= 
BS - 
(=) 
Zz 


Ex 1IBRIS 


Sa Ges 
ra 
Hl __ 
Ss 


DOODSOOSOEEODS 
The University of Alberta 


Printing Department 
Edmonton, Alberta 


re ‘ J 
_ = Cea! 7 cr : 
' ita ay ead) i 
i 7 Gi I M3 


i] : ij 


pg 


7 


ul y en iy 4 rie aete VP 
yi) ies : 
; | eo “ety we dias porns Oy 
ie i | A AD ae Se 
Rebel aeaty a pee titab psc? 


“ We hy a. ¢»" 
/ 7 
i” 


Digitized by the Internet Archive 
in 2023 with funding from 
University of Alberta Library 


https://archive.org/details/Bates1976_0 


THE UNIVERSITY OF ALBERTA 


SLOWNESS-AZIMUTH MEASUREMENTS AND P WAVE 


VELOCITY DISTRIBUTIONS 


by 


ALLAN CLIFFORD BATES 


PO TES ES 
SUBMITTED TO THE FACULTY OF GRADUATE STUDIES AND RESEARCH 
IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE 


OF DOCTOR. OF + PHILOSOPHY 


DEPARTMENT OF PHYSICS 


EDMONTON, ALBERTA 


PALL, W190 76 


DEDICATION 


This thesis is dedicated to my father and mother, 
Matthew and Pearl Bates, my wife, Joanne, and my son, 


Matthew. 


iv 


reson Bas west % S 
,foe wm bas ee ie 2 


We 
(ued 


ABSTRACT 


A study of slowness and azimuth measurements of 
seismic P wave phases, errors associated with these 
measurements, and implications as to the nature of the 
P wave velocity distribution within the earth is 
presented. 

Errors in the arrival times of plane waves at 
seismic arrays result in errors in slowness and azimuth. 
It has been found that the stability of the inversion 
process is related to the condition number of the matrix 
- inversion and that symmetric arrays yield the most stable 
estimates of slowness and azimuth in the presence of 
travel-time errors. If travel-time errors are given in 
the root mean square error sense then the least square 
error. inversion results in maximum and most likely error 
ellipses in the slowness-azimuth error plane for any 
array. The error analysis is not restricted to cases 
for which travel-time errors are specified in the root 
mean Square error manner. In particular if travel-time 
errors are bounded then errors in the slowness-azimuth 
error plane are bounded by multisided figures which 
reflect the array configuration. Berens array confi- 
gurations are used to illustrate that the most effective 
procedure for error improvement is to add additional 
stations along the periphery of an existing array instead 


of in the anterior. 
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Slowness and azimuth observations of teleseismic 
P phases recorded by the 1974 Variable Aperture Seismic 
Array (VASA) indicate departures of the P wave velocity 
distribution within the earth from a spherically symme- 
trical model for which P wave travel times are in 
accordance with the reer are! cuiiey Seismological Tables. 
If lateral inhomogeneities within the crust and upper 
mantle under VASA are not severe, and preliminary indi- 
cations are that this is the case, then the observations 
indicate the existence of anomalously high velocity- 
depth gradients at a depth of about 1200 km, velocity 
inhomogeneities near the core-mantle boundary under some 
areas of the Pacific Ocean, and anomalous conditions 
predominated by lateral velocity gradients between depths 
OL about 1900 and 2600 km under a region close to the 
Caribbean. Also the slowness observations are consis- 
tent with an ‘average’ spherically symmetrical lower 
mantle for which the P wave velocity is slightly greater 
than the Jeffreys-Bullen designation, 

The Tau method of seismic travel-time inversion 
is investigated. A method of estimating extremal values 
bene 0S) paca e 2, 


Re eee a eo where p is the ray 


parameter, T the travel time, and X the epicentral dis- 


of the function 1t(p) 


tance is presented. For each branch of the travel-time 


curve, T observations are fitted to a family of second 


order polynomials in xX. ‘The families of curves are 
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then mapped into the t(p) plane. The Tau method is 
illustrated by inverting T(X) data recorded along the 
"Yukon' line during Project Early Rise. A comparison 
of t(p) results from several other studies in North 
America reflects differences in crustal and upper 
mantle structure. Also the resolving ene a 
function of separation of observation points of the 
Tau method is examined by inverting t(p) envelopes 


calculated from exact velocity-depth functions. 
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CHAPTER 1 


SEISMIC ARRAY DESIGN AND BACKWARD ERROR ANALYSIS 


OF WAVE SLOWNESS AND AZIMUTH 


Introduction 


Array measurements of the slowness and azimuth 
of seismic body wave phases have been used for the 
purpose of delineating earth velocity structure by many 
investigators within the last decade. In this study 
investigations of preers associated with such measure- 
ments and the relationship between the P wave velocity 
of the earth and slowness=-azimuth observations are 
presented. In Chapter 1 errors in slowness and azimuth 
as determined by two dimensional seismic arrays are 
studied using powerful results from matrix iterative 
analysis. The 1974 Variable Aperture Seismic Array 
(VASA) slowness and azimuth observations and the impli- 
cations regarding the nature of the P wave velocity 
distribution within the earth associated with them are 
discussed in Chapter 2. A study of the Tau method of 
travel-time inversion, which is applicable to spheri- 
cally symmetrical velocity functions, is given in Chapter 
3. A method for the determination of errors in the 
function t(p) is presented. The Tau method is illus- 
trated by inverting travel-time data recorded by the 


University of Alberta on Project Early Rise. In all 
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chapters the contributions of past investigations will 
be cited within the text. It is felt that for this 
study such a procedure iS more appropriate than a his- 


torical review. 


The Importance of Array Measurements of Slowness and 


Azimuth 


The two-dimensional seismic array is a major tool 
in the investigation of the velocity and structure of 
the earth's crust and mantle. In most experiments the 
elastic waves are observed in the far field at places 
where the radius of curvature of the wavefront is very 
much greater than the array size and the structural inho- 
mogeneities to be studied. It is possible to assume that 
each phase traverses the array in the form of a plane 
wave with a given azimuthal orientation and slowness. 
If the array is situated on a Cartesian coordinate sys- 
tem, then the equation for the arrival time, t(x,y), of 


a given phase at location (x,y) is, 


thx ey). = eA + Px + PyY 2 Aas 


In (1.1), the kinematic properties of the plane wave are 

specified by tor the arrival time at the origin (x,y) = 

(0,0), and Py and Py? the x and y components of slowness. 
Small seismic arrays are used to advantage in 


exploration geophysics in both the reflection and 
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refraction method to determine mean and pee al veloci- 
ties and to discriminate against noise generated by 
multiple reflections. Medium and large avertnre seismic 
arrays have found wide use in the investigation of ee 
ble departures in upper and lower mantle velocities from 
standard models. A large change in velocity gradient 
may correspond to a small change in velocity and have 
very little effect on the observed travel time. Thus a 
direct measurement of the vector ray parameter Py using 

a seismic array may delineate important anomalies. 
Birtill and Whiteway (1965) have discussed the frequency- 
wavenumber response characteristics of arrays in general. 
Otsuka (1966) gives the expressions for slowness and 
aZimuth variance for given variance of travel-time error. 
The following analysis is a thorough treatment of the 
errors to be expected in determining the velocity and 
azimuth when a plane wave propagates across an array of 
detectors. The analysis is necessary in order to deter- 
mine the array size and shape and the sampling time 
required for a particular significance level. The results 
may be used to determine which observations of anomalies 
in velocity or azimuth are significant and may be inter- 


preted as inhomogeneities within the earth. 
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Backward Error Analysis 


With data from a seismic array, the solution for 


the vector ray parameter, p, iS given by 
paSwoy LT Clo?) 


where A is the ‘array' matrix which depends upon the 
BOGre Males Of thewotations cCompitsingathe array ands 1 
is the 'observation' vector which depends upon the 
observed arrival times and the station coordinates. It 


is common practice to determine the best fitting plane 


wave, in the least squares sense. For N stations situa- 
ted at (X.,Y5), i = 1,N and reporting arrival times tay 
i = 1,N anda best fit plane wave described by equation 


(1.1) we have that the root mean square error, E, in 
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The study of the effects of perturbations of the 
elemertemolethesmacraceA -and-the vector Sl “upon tie eke= 
Nehte=OLevectorcpNishrererred to™in > yumertcal ‘matrix 
analysis as ‘backward error analysis' (Wilkinson, 1965). 
Seismologically this type of analysis is extremely 
important since we would like to know firstly the size 
of the errors in Pp and secondly how to make errors in Pp 
as small as possible. In Appendix 1 some definitions, 
concepts, and results from the numerical analysis of 
linear algebraic systems are introduced (Forsythe and 
Moler, 1967). The results will be applied to the problem 
of measuring the vector ray parameter p. 

In equation (1.2) it is assumed that the elements 
Gepunie «vector Veana possioly the Matrix A are sclrentific 
measurements. For inversion the desired quantity is the 
vector p. Obviously small measurement errors in the 
entries or} fgand T result in errors im /the |}veetorn p. 


it 62 oand BAST represent the errors in T and aT4 respec- 


tively, then we actually measure the quantity p + 6p, 
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It is desirable to design an array in such a way that the 
measured vector p responds ina stable manner, deviating 
from the true p only slightly, upon small perturbations 
in T andy A-ediimrorderjtopobtainisuch ansystemnitiis 
necessary to consider upper bounds of the 'size', in 

the vector norm sense (see Appendix 1), of Op. TE, gas 

is usually the case in seismology, we suppose that there 
aLet aruOons mri (duemt-otervor sttn lannivalytime)sebutyno 


errors in A, then for any vector norm, 
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The ‘condition number', K(A), of a square non-singular 


| . In the 


matrix A is obtained from «(A) = ||A|| [JA 
inequality «k(A) is calculated using the matrix norm (see 
Appendix 1), which is induced by the vector norm used. 
Siti lar bye Cie re: Ales errorse in. bot sAsandi/r then -for 


any vector norm (see Franklin (1968), page 177) 
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The inequalities (1.6) and (1.7) show explicitly the 
importance of the condition number K(A). An optimal 
system results for a given norm when kK(A) is as small as 
possible with respect to that norm; the choice of a par- 


ticular norm depends upon the interpretation and 
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scientific application of the vector p. 


Array Design 


The results (1.6) and (1.7) can now be applied to 
the problem of the determination of slowness and azimuth 
of plane waves. In the slowness-azimuth equation (1.2), 
the elements of the system matrix, A, depend only upon 
the x and y coordinates of the stations. Thus the con- 
dition number, kK(A), depends purely upon the locations 
of the stations. The obvious implication of this fact 
is that we can ensure a stable system, that is minimize 
K(A), by a corresponding correct choice of the geometrical 
pattern formed by the stations. We wish to minimize errors 
in the slowness, Pr and the azimuth, 90. The slowness is 
Saar the Euclidean vector norm, ||p||, of p. The slow- 
ness error is given by [lp+ Sp|| - | |p]. Thus the maximum 
absolute value of the slowness error is ||ép||. From 
(1.6) or (1.7) the maximum slowness error is minimized 
when «(A) calculated using the spectral norm is as small 
as possible. Also, the choice of spectral norm is indeed 
wise for the purpose of designing an array which minimizes 
errors in azimuth. This is true since if 6p is the 


max 


maximum possible slowness error and Cee is the abso- 


lute value of the maximum azimuthal error, then 
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Equation (1.8) is important for two reasons. Firstly it 
shows that a given array is an effective slowness dis- 
criminator if and only if it is an effective azimuth 
discriminator. Secondly equation (1.8) permits a cal- 
culation of maximum ernere in azimuth for given maximum 
error op for any seismic array.. At this point it..is 
worth noting that slowness and azimuth errors, 6p and 
66, cannot achieve their maximum values simultaneously. 
From figure lol it is mot difficult to:see that 

max sp") * 


66 = arctan aes eae e (1.9) 
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From (1.9) we see that 66 attains its maximum value only 
when 6p ~ 0 and that 6p attains its maximum value only 
when 660 = 0. 

Now the inverse problem given by (1.2) will be 
the most. stable when «(A) = 1, where A is given by (1.3). 
By constructing AA and using the results of Appendix 1 
it. is possible to show that «(A) = 1 if and only if the 


following conditions are satisfied; 
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Figure Dee 


The relationship between errors in slowness and 
azimuth is shown here. The vector defined by segment 
OA is the true slowness vector. The radius of the 


error circle is ) oa Three error vectors defined 


Sey 
by segments AD, AC, and AB are shown. Vector AD 


produces the maximum azimuthal deviation OU ie and 


segment AC results in the maximum slowness deviation 


o) SSR The error vector AB contains both slowness and 


azimuth error components; with the length AE=d6p_ the 


relationship between the errors is given by equation 
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The general proof of the above relations is rather 
lengthy. 

The conditions require that the array be cons- 
tructed in a symmetric manner. It is instructive to 
consider the derivation and results when there are N= 3 
stations. “In this case the matrix A in equation (1.3) 


yields, 
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The eigenvalues Sahay of aat are the roots of the charac- 
teristic equation det (aa - oe? = 0orveThus theyssatistfy 
sae IT 7 5 (ath) + (epee } =,O08” Fork (A) Ato besunity, 
the eigenvalues 55: ay must have the same magnitude, and 
hence A), =:11i,“met and sonly ae er = L faebiey4 which 
implies a = b and c = 0. From the above expressions the 


condition c = 0 requires that X= Xo/2. Thius:Gact comp 


bined with the condition that a = b yields Y, Snel: X.. 
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Notice now that the resultant array consists of three 
stations located at the vertices of an equilateral 
triangle. The array locations satisfy equations (1.10) 
and this is the only three station geometry for which 
a: le Pa 

In general arrays which satisfy (1.10) are 
symmetrical with respect to the (X,Y) axes. An example 
of an array which possesses this high degree of symmetry 
is LASA in Montana. Note further that each station of 
LASA is composed of several sensors with the 'sub- 
arrays’ formed by them ‘satisfying (1.10), Thus the 
entire array and each sub-array are capable of reporting 
stable large or small aperture array estimates respec- 
tively of slowness and azimuth. 

It is worth commenting on the effect of array size. 
Suppose we have ne array with system matrix Ay which 
satisfies (1.10). If we increase the size of this array 
by a magnification factor m>1 to form a new array with 


system matrix A then K(A,) =a @urcrons Lor (Chery second 


2° 
array, A,, however, will be smaller since ||A]| and 
[|r] | Wo dbs: dharger,. 

Bor arbitrary Nesp equations. (1. 20))\ tdo not dictate 
a unique geometrical pattern of stations; nevertheless from 


this discussion it is clear’ that it is advantageous to select 


a geometry which does satisfy the above requirements. 
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Error Bounds for Arbitrary Array Configurations and 


Arbitrary Specifications of Arrival-Time Errors 


It is important to realize that the upper bounds 
given by (1.6) and (1.7) in general depend upon the 
azimuth of the incoming plane wave for a given slowness. 
Let us assume that there are no errors in the coor- 
dinates of the stations and focus our attention upon 
the upper bound (1.6). Also suppose for simplicity that 
we have a 3 station determination of p. With the array 
coordinates given by (X.,Y5), i= to) (and for: conven 
ience Xa = ¥, = 0), the upper bound (1.6) can be written 


il 


as follows; 
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The expression (1.11) shows explicitly that the upper 
bound given by (1.6) depends upon the azimuth for an 
arbitrary array. Nevertheless for appropriate choices 
Sfeehesdi rece ons /Ol. the vectors iT and. “oT there’ can be 
equal Cyan 109). ahnus Cor arbitrary. vectors T and .éT 
there is no sharper bound for the relative error 
Wii iP yithaa thet given by (ile), Its for this 
reason that the criterion for array design was based 
upon (16{6)= 

In general the error space defined by vectors op 
for arbitrary arrays.and arbitrary specification of 
travel-time errors ot. is extremely complex. To see 
this consider a three station array with stations at 
(Xi ,¥5), i = 1,2,3, and arrival-time errors of Oty, ét. 
and 6t.,. If we take for simplicity X, = Y, = 0 then we 


5} ai ue 
have. tron CL.3) sand, (1.5.4), 


6p, Y3(6t,-6t,) 2 Y, (6t,-6t,) 
6p -X, (6t,-dt,) + X, (6t,-St,) 


Thus the error space is composed of the space spanned 
by certain linear combinations of three generally 
linearly independent unit vectors Vas ee ee a i 


i and j denote unit vectors along 6p, and SPy axes 


respectively, then the v; are given by, 
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In general then op is a possible error vector if and only 


Zt ele can, be written in the form 
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a3 = 6t,(¥5+X5) “/(X,Y,-X3Y,) : 


Phus rEromrGl aF2 psandel( VN13)k .theherreruspace épyowhich 
depends upon the station locations and the nature of 


the time errors, is given explicitly. Similar expres- 


sions per eaten for (Nex 2) 

Fortunately there is at least one case of travel- 
time error specification for which the error space of 
op vectors is extremely simple even for arbitrary array 
configurations. If there are N stations consider the.N 
dimensional Euclidean vector space with the coordinates 


of the N axes given in terms of the travel-time errors 
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is confined to the hypersphere of 
radius YN 6t, where 6t is a nominal error, then the 
error space in the SE AER 2 plane will be a closed 


ellipse. “le -atis-or tits dase are Given win the Text 


section. 


'Array Error Numbers' and ‘Average' Error 


The diagonal matrix form under orthogonal equi- 
valence (Appendix 2) leads to the simplest interpreta- 
tion about the nature of a square matrix A as repre- 
senting a linear transformation from one Euclidean n 
Space into another such space, It also forms the basis 
bOrpecne: calculations, oL .errors op for arbitrary arrays 
when the: travel—time error vector ot Is confined to the 
hypersphere of radius YN 6t. The resulting error 
equations are extremely powerful for the purposes of 
firstly estimating errors for a given array and also 
for comparing the capacities of different arrays to 
yield stable estimates of slowness and azimuth. 

The results of the preceding paragraph and 
Appendix 2 may now be applied to the determination of 
errors in Pp for any arbitrary array configuration com- 
posed of any number of stations N > 3 for which the 


traveli-timeverro, vector ot = [ory St aeetey is 
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confined to the hypersphere of radius YN §t where 6t 
is a nominal error in arrival time. Extremely useful 
concepts of ‘array error numbers' and average error 
will be introduced. 

The application is possible since we can write 
the error vector ray parameter equation for N > 3 in 


the following form; 


Me aaG a ti (ivi) 


where SP aug is an ‘augmented! (1 x N) error vector ray 


parameter with (N-2) zeroes given by, 


iE 


8Paug = (6PyroPyr0r0,---,0) 


and B is an ‘augmented!’ array (N x N) matrix which has 
all zero entries in its (N-2)-th bottom most rows. The 
non-zero top two rows of B depend only upon the station 
coordinates. The matrix B has two and only two non- 
zero Singular values Hy? Uo hence B maps hyperspheres 
in _ §t.space.into ,ewo, dimensional: ellipses; in Paug 
space. The diagonal orthogonal equivalent of B, namely 
D say, is of Ene, form. Dw eed ce) Itievewilhainene 


1) 
Li aanpin’ Hea pls and all other a5 = 0. The orthogonal 
transformation V (Appendix 2) merely 'rotates' the 6ét 
space leaving any hypersphere of radius YN 6t containing 


the same elements as before. Also, it is not difficult 


to show that the orthogonal matrix U (Appendix 2) 


aly) 


“otc Rak ee ody ~ : 
“mt ee tra nl 


or a Ses 


»sa~ 


Oro. » . . soa = i a 
salt bh my A 
VRE 3edomw sox19 u « £) "Bed nompts' rnss al usd? © 


Yd Weavip Bs (Si) coiw 


ee qd) = eae nd 


rae 


eed dokdw xitgzem (Wx WW) a aga ‘besaempus' 16 wh a 
~ De 


sat ,awox yaom mottod i ais eat mk eaeixtne aa £. 


wokiste ons croqus yino . Grahactal 4 to awox ow? qos © Su 
~ton owt yino a ows esr a kiatsm edT 


$3 
wail? mt edegifis pidtaie dais ows odnmi oosae “2m 


ae ,a 20 tne levivpe fsanoportro Lsnopsib of? 
oneriw (Wy..aygl = Ege we = 4 miot ofa Io ek 
{snopons1o ont 0 = eed saris fis bag ou =, gel “so 4 ae 
| 22 of ‘sogetox' yietom is A) V noes: ie Be 
meeieepte a 
2a ee eee -Suoted a6, atnemele 


‘ eibeon 9 nse Lenoanao ot Ja 


a 


operates on the Soe kig space in such a manner that it 
only rotates the space in the (Soar plane. Thus, 
if we make an observation Detox which the root mean 
square error in arrival time is §ét then the error 
vector op will be contained within an ellipse in 
(Sp, Sp ) space which has semi-major and semi-minor 


y 
_ axes, Op and OP respectively given by, 


6p) = H, YN 6t 


Cire BS) 


6P5 = Uy NOE OM 


The quantities 6Py and Op, may aptly be referred to as 
rarray error numbers’. They represent maximum possible 
6rrors for a given ‘root mean square’ error, st. in 
arrival time. Among the vectors 6t of length YN 6t 

Meo Gesvacc there wills be oneadirearionsror which ot 
Bude octemap Onto thewtwo exuremlt1es Of theysemi—major 
axes of the sopyerror ellipse; there will be a direction 
orthogonaly to Chis one for which ot and 0c map onto 
the extremities of the semi-minor axes of the ellipse. 


The other (N-2) mutually orthogonal directions are 


mapped into Sp = 0 for any length ||ét||. The quanti- 
ties Spy and 6p5 are useful since they do give the 
maximum errors, that is the extremities of the ellipse. 


For the purpose of comparing errors admitted by two 


different arrays with the same number of stations one 
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could compare the areas, TOP, OPor of the associated 
maximum error ellipses of the two arrays. 
Points on the maximum error ellipse described 


above are damage points: ofivectorss 6t. of: length YN 6t 


which lie entirely within one distinct hyperplane in the 


6t error space. As the number of stations increase.the 
dimensioniof, the Dt error.wspaceéiuncreasess andthence* the 
preobabilityithat Ot ties entirely within any one hyper- 
plane decreases. Thus the 'effective' errors not only 
depend upon the maximum errors OP, and OP but also 

upon the number of stations. In order to study the 

effect, ofthe mumber)j,of) stationselétp thei stterrorpspace 


be orientated such that it is spanned by unit vectors 


ot,, 1 = 1,...,N and 
T 
B(YN dt 6t,) = (dp,,0,...,0) 
B Wi fits dn 200 SP5 Oopacs 0)” 
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SOP aLMLamVecLOr sot. nel Lemgi no F Op = Bdédt depends 

Upon EL ne Component (Of sotwi ln. the (ot, , St.) plane which 

is (85485) %, Peesthe none ae Seine Sey a NG errs 
equally likely random variables subject to the cons- 
cratntethac ZB = Net, then the most likely value of 

gf + Bs will be WS Also, Since the direction of the 
Projection Of soln the (ot, , St.) plane is random, the 
Slireccave Oreaverage Grromicpacesor vectors (6p will ‘be 
defined by the image of a circle in the (St, ,5t,) plane 

of radius (87485) = V/27 ote. under the transformation B. 
This image will simply be an ellipse in the (Sp. 9P,;) 
plane with semi-major and semi-minor axes lengths given by 
(2/N)* 6p, and (2/N) * 8p, respectively. Notice that this 
‘average’ ellipse is simply a reduced version of the 
ellipse describing maximum errors. A realistic quantity 
which reflects an array's capacity to discriminate against 
random errors is certainly the area of the above ‘average! 
-ellipse which is’ (2176p, 6p,)7N. ~Alternatively~ one-could 
consider the radius of the circle which has the same 


area as the given ellipse, namely, 
1 
Ue 3 
6p, = ((26p, 6p5)/N}* . Cr 16) 


The value OPy can be referred to as an ‘average error' 
Since it takes into account the number of stations and 
the random nature of the true errors and it is an average 


over the skewness (elliptical shape) of the error space. 
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The values 6P119Por and SPy depend only upon the 
station coordinates. The singular values Uy and U5 of 
B are the positive square roots of the two (necessarily 
positive) non-zero eigenvalues of the matrix C = ppl, 
TaegeneraisiCa will be of sthe.torm (C i= less, deg 0 eee eee NL 
where Oa OU foieo> GSVOm fo jon. nengenerally non-zero 
depend only upon the 


quantities c and c 


Bal Wine 2:4 ie Bon 
station coordinates. The quantities Uy and U5 are given 
by; 
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* 2 %% 2 
Hy = {ley yteggtl (cy +Cg9)° - 4,4 Cn9-Cp9)1 1/23 


Cha) 


i 2 2 % 
Hy = Eley tegg-[ ley +Cg9) - 4 (C44 So9-C7 9) 1)/ 23. 


POMse Oya bop and. (is) othe anray. error numbers 6p, 
and OD5 are given. The average error OPa is then given 
Dyeeh. 16% 

hoy LoenOt. Glerrowvemcto, mind, the “Eile of tne 
error ellipse in the Greer ee plane (Golub and Kahan 
(1964)). Recalling the results of Appendix 2 we have 
that (8D, Py 1042-20)” = U(6py, 6B) ,0,.--,0)°. Now the 
columns of the orthogonal matrix U are simply the eigen- 
vectors of the matrix C; this fact permits the calcula- 
tion of U and hence the angle of rotation in the 
(Spr Spy) plane. Note also that Golub and Kahan (1964) 


have devised a computational procedure for determining 


the singular values of general matrices; the procedure, 


7 £ 4 7 
Lie SE 


_ é. + 
: : o.oo a oo oe ‘ 
isang YO esutny el pe 
i eae? y bas) « 
te 


>. Pe ‘ 
— . 
_ eu 


Tae 
+ 


which is extremely complex analytically, was introduced 
since in general the calculation of BB- for arbitrary 
matrices B, using floating point arithmetic does serious 
violence to the smaller singular values as well as the 
corresponding eigenvectors which appear in the matrices 
U and V of Appendix 2. It is felt, however, that for 
this analysis the treatment of Golub and Kahan (1964) 
is not necessary since the matrices involved are ex- 
tremely simple (for example, B contains only two non-zero 
rows). 

Now we need to determine expressions for the 
augmented array matrix B, the quantities OP, 1 9P5 and 
hence OPa- Manipulation of equations (1.3) yield B. 
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The values 6P,19P5 and Opa are given by (1.15), (1.16) 
and+1 7) -""1t"is* important *to irealizge-that ‘these 
quantities do not depend upon the location in the (X,Y) 
plane from which the X; and Ys are measured. 

In a previous section it was stressed that the 


most stable estimates of p result when equations (1.10) 
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are, Satistiied. It is not GLfticult to show that if 
(1.10) are satisfied then the 6p error space becomes 
arcarcle: with 6p) = 6P5 and SPy = (2/N) * 6p, 

It has been shown that if travel-time errors 
combine to give a root mean square error of YN 6ét 
then the resulting maximum errors in p will be such 
that 6p has its tip anywhere on the circumference of 
an ellipse centered at the origin in the Op oP) 
plane. Also there will be another ellipse, with semi- 
major and semi-minor axes lengths dependent upon the number 
of stations, which defines 'most likely' errors Op. 
Examples of the maximum and most likely error ellipses 
for several array configurations will be given in the 
next section. It is also interesting to consider 
resultant maximum errors in the (SP. OP) plane when 
travel-time errors are not described in the root mean 
square error sense but in the 'norm-infinity' sense. 
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NXN matrix A, the matrix norm infinity is defined by 
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| |A| |. = max y |a, .| where the a, . are the 
fees Niwa dete ae 
entries of the matrix. Thus if ||ét||, is given one 


could determine the quantity ||ép||, = max{| dp, |,| 6p, |} 
using |jiep| |, <= | (B)IJtiot}|. and equations (1.18). 
Instead of considering this maximum, ||ép||_, it is much 


more interesting and instructive to determine the 
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"boundary' in the (Spy + SP,,) plane when it is assumed 
that -s6t < ot; = Ot: for, aly ét 5 where §t > 0 iS a 
nominal travel-time error. Notice that this travel- 
time error specification is quite different than the 
root mean square error case. Using equations (1.14) 
and (1.18) it is possible to determine this boundary. 
In general the boundary in the SPL 7 SP,,) plane is a 
multisided figure which has a degree of symmetry 
dependent upon the degree of symmetry of the array. 
This maximum 'norm-infinity' Pome will also be 


shown in the following examples. 


Examples 


Much emphasis has been placed upon the concepts 
of condition number, maximum and most likely error 
ellipses, and the average error OP. The following 
examples illustrate the importance of these concepts 
and also show the relationship between spectral or root 
mean square errors and norm-infinity errors... In order 
to compare the capacities of different arrays to dis- 
criminate against travel-time errors it is necessary 
to ensure that the sizes, in some sense, of the arrays 
are comparable. To this end the 'mean' aperture of all 
of the following arrays are equal. The horizontal 
extent of an array depends upon the azimuth from which 


the array is viewed; if ay and a, are the largest and 
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smallest apertures of the array then the geometric mean 
aperture, a_, TondeLineditompe alties Va a.- Figures 
Dee ee esp alt wane). 5) SHOW various array configurations 
together with their corresponding spectral and norm- 
infinity error spaces. The distance scale factor AL 
and slowness scale factor 6p are such that when AL = 
200 km and the nominal travel-time error 6t = 0.1 sec, 
then 6p = 0.1 sec/100 km and all of the arrays have a 
mean aperture of 200 km. The errors for any travel-time 
Berean array size can thus be determined by using the 
relation ép =< 6t/AL. 

Three, five, and nine station L-shaped arrays 
albenonown win Lugures 2a.) LJ2b, and 1 o2c; they twill 
be referred to as arrays 2a, 2b, and 2c. The feature 
that all three of these arrays have in common is that 
they are asymmetric; accordingly for all three arrays 
the spectral error space is elliptical, not circular, 
and K(A) > 1. The tilt of the ellipses (the semi-major 
axis is inclined at an angle of 45° counterclockwise 
from the op, axis) reflects the asymmetry of the arrays. 
Notice that as we proceed from three to five to nine 


stations the following take place; 


(i) Enieweondst10on numbers ancrease (from 2.5 to, 2.9 
Oe 3.1) 
(4-7) the spectral error spaces become more elliptical 


(error ellipses acquire larger eccentricities). 
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Figure dee 


Three, five, and nine station L-shaped arrays are 
shown in parts a, b, and c respectively. Immediately 
above each array the maximum spectral error ellipse 
(light boundary) and most likely spectral error ellipse 
(dark boundary) are shown. The scale factors shown are 
such that when L = 100 km, and ét = 0.2 sec then 6p = 
AO oe sec/km. In this case the values of 6p, (in 
the order a, b, and c) are 6p, = 4.6 x 100 sec/km, 


3 3 


4.3 x10 ~ sec/km, and 3.6x10~ sec/km. If L = 100 m 


and é6t = 2 ms, then dp = 4.0 inna sec/km and (in the 


2 iz 


order a, b, c) 6p, becomes 4.6 x 10 * Ssec/km, 4225-7 200 
sec/km, and 3.6 x 100- sec/km. Immediately above the 
spectral error space diagrams the maximum norm-infinity 
error boundary is shown for all three arrays. The hash 


marks on the slowness axes of all slowness error diagrams 


are Op from the origin. 
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Figure i103 


Three, five, seven, and nine station "circular" 
arrays are shown in parts a, b, c, and d respectively. 
The corresponding maximum spectral error circle (light 
boundary) and most likely spectral error circle (dark 
boundary) are shown immediately above each array. The 
scale factors shown are such that when L = 100 km, and 
6t = 0.2 sec, then ép = 4.0 ale sec/km. In this case 
the values of ODay (in the order a, b, c, and d) are 
4.6x10°° sec/km, 4.2 x107° sec/km, 3.8 x10"> sec/km, 
and 3.4% 1057 sec/km. If L = 100 m and 6t = 2 ms, then 
6p = 4.0*~x ne sec/km and (in the order a, b, c) 6Py is 


AvGou10;- secVkm, 402 2108 ceca Seon anome 


2 


sec/km, and 
3.4 x10 “ sec/km. Immediately above the spectral error 
diagrams the maximum norm-infinity error space is shown 


for all four arrays. The hash marks on all slowness 


error diagrams are at 6p from the origin. 
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Figure 1.4 

Two nine station and a thirteen station "double 
circle" arrays appear in parts a, b, and c respectively. 
The corresponding maximum spectral error circle (light 
boundary) and most likely spectral error circle (dark 
boundary) are shown immediately above each array. The 
scale factors shown are such that when L = 100 km and 
ét = 0.2 sec, then Sp = 4.0 VAs sec/km. In this case 
the values of SP, (in the order a, b, and c) are Opp = 


3 3 sec/km, and 3.2 x107> sec/km. 


37 10%. esec7km;, 32/0 
TfL = 100em and st = 24ms, tren ope= 4. 0r 10° sec/km 
and (in the order a, b, c) 6Py LS. 37 Sane sec/km, 


# sec/km, and 3.2 108 sec/km. Immediately 


207 Xow 
above the spectral error diagrams the maximum norm- 
infinity error space is shown for each array. The hash 


marks on all slowness error diagrams are at 6p from the 


origin. 
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Figure 1.5 


Nine, seventeen, and twenty-five station "box" 
arrayS appear in parts a, b, and c respectively. The 
corresponding maximum spectral error circle (light 
boundary) and most likely spectral error circle (dark 
boundary) are shown immediately above each array. The 
scale factors shown are such that when L = 100 km, and 


3 


6t = 0.2 sec, then ép = 4.0 x10 ~ sec/km. In this case 


the values of SPy (in the order a, b, and c) are 6Pa = 


3 3 


3.4*x 10 ~ sec/km, 2.5 x10 sec/km, and 2.4 x10 ~ sec/km. 


If L = 100 m and ét = 2 ms, then 6p = 4.0 x107* sec/km 


- sec/km, 


and (in the order a, b, c) dp, is 3.4 x10 
Py So aay - sec/km, and 2.4 Om sec/km. Immediately above 
the spectral error diagrams the maximum norm-infinity 


error space is shown for each array. The hash marks on 


all slowness error diagrams are at 6p from the origin. 
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(iii) the arrays become more and more asymmetric with 
the L-shape being more well defined. 
The above features are related. In general larger con- 
dition numbers are associated with more asymmetric 
arrays and more eccentric error spaces. One might 
suspect that the ratio a./ay is equal to the ratio of 
the semi-minor to the semi-major axes of the ellipse. 
This is not the case however; for example using array 
2a we get a./a, = 0.5 and the ratio of the ellipse axes 
is approximately 0.57. The reason for this is that the 
analysis not only takes into account the extreme 
apertures of an array but also the spatial distribution 
of stations along the directions of the extreme aper- 
tures. Notice also that as the number of stations is 
increased the maximum error ellipse becomes larger. 
This is a result of the fact that the addition of 
Sone is performed in a manner such that the station 
spacing is reduced. On the other hand, the most likely 
error ellipse diminishes in extent as the number of 
stations is increased. Starting with array 2a the 
improvement obtained by the addition of 2 stations 
(array 2b) is such that Spy decreases by approximately 7%. 
Four additional stations (2c) improve errors by 16% when 
compared to array 2b. Also shown is the maximum norm- 
infinity error space; the boundary is generally multi- 
sided and its orientation again reflects the asymmetry 


of the array. As the number of stations increase the 
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boundary acquires more facets approaching the elliptical 
shape of the spectral error space. The boundary of the 
es Oe error space falls between the maximum and 
most likely spectral error ellipses in all cases. 

Pigure b.3 shows! 3, 5; ./ and.9 ‘station cireular 
arrays and their error spaces. All of these arrays have 
station locations which satisfy equation (1.10) and 
hence kK(A) = 1 and the spectral error spaces are cir- 
cular. The radius of the maximum error circle increases 
as the number of stations increase. The radius of the 
most likely error circle, as expected, decreases as the 
number of stations increase. In fact Opa decreases 
almost linearly as the number of stations increase; ODay 
diminishes by approximately 6é6t/5AL with the addition of 
every two Stations. The norm-infinity error spaces are 
multisided figures with a higher degree of symmetry than 
those for the L-shaped arrays. The norm-infinity boundary 
for the nine station array. 3d is the most !circular'. 
Again the norm-infinity error boundary lies between the 
maximum and most likely spectral error circles, It is 
interesting to compare the effectiveness of these circular 
arrays and the previously discussed L-shaped arrays. 
Comparing arrays with the same number of stations it is 
seen that the circular arrays are more favourable in all 
cases with the improvement in corresponding SPy being 


the largest for the nine station arrays. The disparity 
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in SP is not extreme with the ratio of 6Py for array 

3d to SPy for array 2c being approximately 0.94. A 
comparison of maximum errors allowed by the two array 
configurations again shows that the circular arrays are 
more favourable in all cases and the improvement is 
again most pronounced for the nine station arrays. The 
disparity in maximum errors is much more pronounced than 
the improvement in Opai for example the ratio of the 
maximum error allowed by array 3d to that of array 2c 

is approximately 0.67. The superiority of the circular 


arrays is related to the fact that their corresponding 


condition numbers are equal to unity; the L-shaped arrays 


have K(A) > 1. 

Figure 1.4 shows the double circular arrays 4a, 
4b, and 4c. Notice that the spectral error spaces of 
the nine station arrays 4a and 4b are identical. This 
‘rotational symmetry! is a feature common to arrays for 
which «k(A) = 1. Notice however that the maximum norm- 
infinity errors of 4a and 4b do not coincide. A com- 
parison of the value OPy for arrays 4a and 4b with the 
value of Opp for the nine station array 3d shows that 
the single circle arrangement is preferable; OPy is 
smaller by about 10 percent. The thirteen station array 
4c shown has a configuration which contains many equi- 
lateral triangle components which are themselves capable 


of smaller aperture estimates of slowness and azimuth. 
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Its maximum norm-infinity boundary is multisided and 
closely approximates a circle, 
| The nine, seventeen, and twenty-five station box 
arrays are shown in figures 1.5a, 1.5b, and 1.5c res- 
pectively; all have circular spectral error spaces. 
The maximum norm-infinity error space becomes more 
"circular' and larger as the number of stations is 
increased. A comparison of the effectiveness of these 
arrays is very interesting. The addition of 8 extra 
stations to the 'periphery' of array 5a results in array 
5b; the error improvement is such that OPy decreases by 
about 26 percent. If an additional 8 stations are added, 
this time to the interior of the box, the result is array 
5c. This last step results in an improvement of only 
about 7 percent. This situation again shows that an 
array is most effectively upgraded by the addition of 
stations if the new stations are placed along the peri- 
phery of the existing array and not added to the interior. 
A general feature of spectral errors is that the 
size of spectral error spaces is directly proportional 
to the nominal travel-time error 6t and inversely pro- 
portional to array size. The travel-time error depen- 
dence is clearly seen and the nature of the array size 
dependence can be explicitly stated as follows; if an 
array defined by (Xi ,Ys), i ie wee Ne DAS. errors 6P1719Po 


and 6Py then the array defined by (cX,,cY,), cal Ama rs | 
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where c>0O has errors P)/c, Po/c, and Pa/C- Now 
let us consider a 21 element L-shaped array for which 
each leg of the L is 25 km and accordingly the station 
spacing along each leg is 2.5 km. If 6t = 0.1 Sec, 
then 6p, = 2.3 x107* sec/km, 6p, = 1.2 x10°* sec/km, 


and ép, = 0.5 x 1072 


sec/km. If this array is now 
magnified by a factor c = 8 such that the legs of the L 
become 200 km and the station spacing along each leg is 
20 km, then the errors are (for é6t = 0.1 sec) 6p, = 


3 3 


2.9x 10 ~ sec/km, 6po = 1.4 x10 ~ sec/km, and Spay = 


0.6x 107° 


sec/km. Now an important travel-time error 
results from miscorrelation of phases from one record 
to another; in general the accuracy of correlation 
increases as station spacing becomes smaller. Thus if 
phase identification is the only source of error, then 
the smaller L-shaped array described above will be as 
effective as the larger array if the smaller station 


spacing decreases the phase identification error by a 


factor of 8. 


Conclusions 


Errors in slowness and azimuth measurements from 
two dimensional seismic arrays depend upon array geome- 
try, array size, and the nature of travel-time errors. 


For travel-time errors defined in the root mean square 
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error sense it has been shown that the maximum errors 

in slowness and azimuth are defined by a spectral error 
ellipse in the slowness-azimuth error diagram for any 
array. Values of expected errors are then defined by 

a reduced vereiae of the maximum spectral error ellipse. 
In general the stability of the inverse problem is 
related to the condition number associated with an array. 
Symmetric arrays have associated condition numbers equal 
to unity and circular spectral error spaces in the 
slowness-azimuth error plane. AS array asymmetry becomes 
more pronounced, condition numbers increase and spectral 
error spaces become more elliptical. Comparisons between 
symmetric and asymmetric arrays with the same number of 
stations and equal mean apertures reveal that symmetric 
arrays are more favourable in terms of stable measure- 
ments of slowness and azimuth. The equations developed 
show that slowness and azimuth error analysis need not 

be restricted to the case of spectral errors. In par- 
ticular if it assumed that travel-time errors at each 
station are bounded, then it is possible to construct 

the 'norm-infinity' error boundary in the slowness- 
azimuth error plane. The norm-infinity error space 
boundaries are generally multisided figures with a degree 
of symmetry dependent upon the degree of symmetry of the 


array and the travel-time error bounds. 
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The analysis shows that the most effective proce- 
dure for error improvement is to add additional stations 
along the periphery of an existing array and to make the 
array aS symmetric as possible. The size and station 
Spacing of the array should, of course, be as large as 
possible but must be limited by the spectral character 
of the noise, the degree of correlation between signals 
at widely spaced stations, and the validity of the plane 
wave approximation. A large number of stations may be 
necessary for enhancement of weak signals in the presence 
of strong random noise but if the signal to noise ratio 
is adequate and the errors are due to terrestrial inho- 
mogeneities it is probably better to use several adjacent 
arrayS with a modest number of detectors (5 to 9) instead 


of a Single array with many detectors. 
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CHAPTER 2 


TELESEISMIC SLOWNESS AND AZIMUTH MEASUREMENTS : 


THE 1974 VARIABLE APERTURE SEISMIC ARRAY 


Introduction 


Medium and large aperture seismic arrays have been 

idely used in past for the determination of the apparent 
slowness and azimuth of various teleseismic phases. In 
this chapter a study of the slowness and Apia here? 
phases, for which the rays bottom in the lower mantle, 
as recorded by the 1974 Variable Aperture Seismic Array, 

(VASA), will be presented. Usually results are compared 
to values predicted by some standard earth model such as 
that given by Jeffreys and Bullen. The Jette, s-Hul len 
velocity depth profile is spherically symmetric and is 
consistent with an earth that.is chemically homogeneous 
within the lower mantle. A seeaeture of the velocity 
gradient at a given depth from ‘ena given by the 
Jeffreys-Bullen model points to an ‘inhomogeneity! in 
the lower mantle. Changes in the ray parameter, p, as 
a function of distance are extremely sensitive to 
anomalous velocity gradients. On the other hand, a 
large change in velocity gradient which corresponds to 
a small change in velocity has very little effect on 
the observed travel time between about 30 and 100 degrees 


epicentral distance. Thus array measurements of p= daT/dA 
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for the P phase are extremely important. In early studies 
it was assumed that departures of the velocities within 
the lower mantle from the Jeffreys-Bullen model were 
radial in nature. Thus, initially it was assumed that 
the azimuth of the incoming plane wave was the great 
circle azimuth between the array and the location of the 
seismic event as given by the world-wide seismic network 
determination; calculations of p were made with this 
azimuthal constraint. In later studies Gopalakrishnan 
(1969) for example, calculations of both slowness and 
azimuth were made but only information about the ray 
parameter p wasS retained; any event for which the cal- 
culated azimuth deviated from the expected great circle 
azimuth by more than about 4 degrees was discarded on 
the basis that the measurement was grossly in error. 
Still later it was thought that the lower mantle could 
possibly be laterally as well as radially homogeneous. 
Thus both slowness and azimuth information were retained. 
Davies and Sheppard (1972), for example, plotted slowness 
and azimuthal deviations on the ‘array diagram’. Thus 
both slowness and azimuth measurements are important. 

In accordance with this, results from VASA 1974 will 

be presented in the standard A(p) plane and in the slow- 
ness-azimuth ‘array diagram’. Also in view of the 
importance of the function t(p) (Bessonova et al (1974) 


and Bessonova et al (1976)) in crustal, upper mantle, and 


eins itiw oben oew a to, snolseinals® tol Jamie Ies 
ignidat satel egoo eoibuda 3038 os axtanoo f , 


Vern) ots auc + naka 4 i 

~{Leo ois doidw sot 109V8, OR tis 

slowto tsexp, bstoegke ant moh. | 
ng bebusneiS asw aaexpsb' bh porns asta om nba MLSS 

oh ¥ aie ; 

.touxs at yleeorp 2Bw themexwasom ons vais Waal a : . 

xives olinsm xewol ent Jans stgvans aa 9 Raainy Fr oN 

-eposmmpomod yt Lakpas 26 [isw #6 viseresel on 


.bethetet exew nok teirsoha djumtee brs seonwole titod | 
A Saws 


seonwelea bottolg siqmpxe ae tenes) Brsccerta - 
sud? .'msipeib yerzs’ ext ite) enoksobveh tnd 
-tas2t3toqnl 926 adnomazvesen 


iftw BvCL ABAV mox? eviuees aye 


x hin 


-wole ads ai bos oredg rT Po eS 
ans 30 weiv ni oatn el pan 


(NGL) Ee te evonoaee€) (qd £soe ‘ 


ew 
re she ; 
or mi” ® 


44 


lower mantle studies, results will also be presented in 
the t(p) plane. Before these results are presented, 
details of the 1974 variable aperture seismic array, the 
method of data acquisition and reduction, and the 
"COVESPA' process (used to determine slowness and azimuth) 


will be given. 


The 1974 Variable Aperture Seismic Array 

The portable 1974 VASA consists of the six most 
southerly stations shown in figure 2.1. Notice that this 
arrangement of stations is highly symmetrical, an impor- 
tant feature in view of the results of Chapter 1. In 
order to increase its capacity, the array was extended 
to include the permanent station EDM. For the purpose 
of calculating dT/dA aes azimuth of teleseismic P waves, 
this extended array was divided into two sub-arrays, 
VASA1 and VASA2. The sub-array, VASAl1, consists of the 
stations EAT, FOR, HAN, MAP, SES, and VUL, and VASA2 
consists of the stations EAT, HAN, SES, VUL, and EDM. 
Thus the arrays VASAl1 and VASA2 are large aperture type 
seismic arrays and are therefore capable of yielding 
stable estimates of dT/dA and azimuth. For reference, 
the locations and elevations of the stations are given 


in Table 2.1. 
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Figure 2.1 


Array location map. The 1974 Variable Aperture 
Seismic Array, VASA, consists of the six most southerly 
stations Eatonia (EAT), Foremost (FOR), Hanna (HAN), 
Mape Creek (MAP), Suffield (SES), and Vulcan (VUL). 

The permanent station Edmonton (EDM) augments the 1974 


VASA. 
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Table 2.1 


1974 VASA Station Locations 


Code ' Elevation 
Station Name Latitude Longtitude (km) 
Eatonia, EAT De OO LS O3 O2:70 


Saskatchewan 


Foremost, Alberta FOR ASO VS5as40 21274430733" 0.97 


Hanna, Alberta HAN Sree 7 205) tee” 3.73" 0.93 
Maple Creek, MAP 49° 47.73% 109° 20.52' 0.95 
Saskatchewan 


Suffield, Alberta SES SOP 2s.) ete’ s 2750" O01 7 
Vulcan, Alberta VUL S02 e505 Leslee cic oe L007 


Edmonton, Alberta EDM 53° 13.37' 113° 20.90' 0.73 


Each station consists of three Willmore Mark II 
seismometers (one vertical and two horizontal measuring 
north-south and east-west motion), a WWVB receiver, 
amplifiers, a multiplexer, and an analog to digital 
converter. Kanasewich et al (1974) discuss the essen- 
tial features of the tripartite digital recording gain 
ranging system. The response of the amplifier and a 
Willmore Mark II seismometer is shown in figure 2.2. 
The desirable feature of the response is that it is 
flat over the frequency range of teleseismic P waves 
(about 0.7 to 1.8 Hz). 

The procedure for editing raw seismic data recorded 


by the portable field system has become routine at the 
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Figure 2.2 
Amplifier gain of the tripartite system and 
the combined response of the amplifier and a Wilmore 


Mark II seismometer. 
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University of Alberta. Data selected from the 7-track 
magnetic field tapes are transferred to 9-track tape 

by means of a PDP 11 computer. In general the sampling 
rate is variable but data used in this experiment were 
sampled at 12.5 samples/sec per channel (except for a 
few at SES which were sampled at 65 samples/sec per 
channel). Information from each seismic event, for each 
station, is contained in two files on final master tapes. 
The first file contains the statistics of the event; 
that is the location of the epicenter, the magnitude, 
the focal depth, the origin time, the distance and 
azimuth from the station, the Jeffreys-Bullen arrival 
time of the P phase at the station, and the ellipticity 
and elevation corrections. The second file contains the 
seismic data - the vertical, radial, and transverse 
traces and the WWVB signal (radial and transverse traces 
are obtained assuming a great circle path from the epi- 
center to the station). All four channels are stored in 
(four byte) integer format with the data file composed 
of blocks of 8192 bytes. With a sampling rate of 12.5 
samples/sec per channel this means that there are 
approximately 163.84 seconds of seismic information 

per block. The data on the final master tapes is 
written in cyclical fashion with respect to the channels; 
thus the first four words are the first entries of the 


vertical, radial and transverse traces and the WWVB 
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Signal respectively, with the fifth word being the 
second entry of the vertical channel. For a sampling 
rate of 12.5 samples/sec per channel, there is a time 
shift between successive channels of 20 milliseconds 

on the original field tapes. On the final master tapes 
the first three channels are time shifted, by means of 

a simple linear interpolation scheme, so that all three 
are synchronous with the WWVB channel. Accurate time 
resolution is accomplished by cross correlating the 
actual WWVB signal with a synthetic WWVB time signal in 
binary pulse code format 96 sec long. A visual deter- 
mMination of the actual WWVB signal ensures that both 
Signals contain the same minute mark. The maximum 
cross-correlation then corresponds to the coincidence 

of the two minute marks and the fraction of the sampling 
interval (80 ms) for which the maximum occurs is deter- 
mined by a three point fit to the correlation values. 

It has been pointed out by Gutowski (1974) that for good 
Signal to noise characteristics the precision of such a 
time determination is to within 1/4 of a sampling inter- 
val which is 20 ms for a sampling rate of 12.5 samples/sec 
per channel. Verification of this result has been per- 
formed by decoding separated portions of the actual WWVB 
Signal. With accurate timing available the master tapes 
have been arranged such that the time corresponding to 


the beginning of the first block is two minutes before 
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the JB arrival time of the P phase. Data on the master 
tapes are written with all stations recording an event 
following one another. Events are arranged chronolo- 
gically according to their time of occurrence. For 

this experiment the data from the station EDM required 
special handling. The EDM data is transmitted via a 
data link to the seismic laboratory at the University 

of Alberta with End rate of 18 samples/sec per 
channel. The EDM data has been modified so that it is 
present on the master tapes in the standard form. As 

in past (Gutowski (1974)) the organizational scheme 
described above has proved to be convenient for the study 
of teleseismic P phases. Also, it is expected that such 
an arrangement will prove useful for future seismic 


studies. 


Velocity Spectral Analysis and the Covespa Process 


Velocity filtering techniques have been widely 
used for the determination of the apparent velocities 
of wavefronts which traverse seismic arrays. The design 
of many velocity filters is based upon the assumptions 
that the frequency spectra of signal and noise are dis- 
joint or that the noise, unlike the signal, is incoherent 
across the array. An exhaustive treatment of the time 
delay - summation process (stacking) and cross correla- 


tion techniques, in conjunction with one and two 
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dimensional arrays, is given by Birtill and Whiteway 
(1965). In the oil industry, the examination of contour 
plots of cross correlation values or 'coherency' in the 
apparent velocity-time plane has long been standard 
practise for the purposes of the identification of 
primary and multiple reflections and the determination 

of interval ‘LER ME (Gee and Backus (1968) and 
Taner and Koehler (1969), for example). Davi ee et al 
(1971) used the summation process to determine the 
velocity spectra, 'Vespa', of teleseismic wavefronts 
traversing LASA; the azimuth of the wavefronts was 
assumed to be that given by the great circle path between 
the event and the array and results were illustrated by 
the 'Vespagram'. Later, Doornbos and Husebye (1972) used 
the method of Davies et al (1971) and cross correlation 
techniques in order to study the apparent velocity of 
core phases. Wiechert et al (1967) determined the velos 
city and azimuth of teleseismic waves arriving at the 
Yellowknife seismic array by examining the squared summa- 
tion of traces. A complete treatment of velocity filter- 
ing in general is given by Kanasewich (1975). 

The 'Covespa' technique (Gutowski (1974)) is an 
extremely powerful spectral method ops determining slow- 
ness and azimuth and has been employed in this study. 

For a given azimuth, 0, slowness, p, and time, t, the 


coherency, CC, is given by, 
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2.1) 
i 


2 asses 
CC(6,p,t) = ue ) ) ) t 
M(M-1)-T ere {} ¢2 


where M is the number of sensors, k is an incremental 
integer on channel i (i # k), T is the length of the 
time window (1 second for this study), and f. is the 


,c 
th channel at time t. The traces from 


amplitude of the i 
the M sensors are shifted in time by the amounts deter- 
mined by the particular choices of p and 6. For each 
time along the records the zero lag cross correlations 
of all combinations of two stations are computed, nor- 
malized to unity and summed. Thus CC = 1 at a given 
time, Slowness, and azimuth if the phases and waveforms 
of the signal within the window are the same at all 
sensors. For this study the Covespa technique has been 
applied to the vertical traces with a time window of 1 
second. 

A contour plot of coherency in the slowness-time 
plane is known as a 'Covespagram'. The range of accep- 
table coherencies is set at 0.5 < CC < 1.0 to ensure 
meaningful interpretation of Covespagrams. Gutowski 
(1974) has shown that the Covespagram pattern remains 
stable upon variation of relative amplitudes of wave- 
forms across an array and that an increase in wave train 


duration merely results in an extended Covespa pattern. 
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An example of the resolving capacity of the Covespa 
process when applied to coherent signals embedded in 
"random' noise is here shown in figures 2.3 and 2.4 

for signal to noise amplitude ratios of 4:1 and 1:1 
respectively. The signals are six synthetic two cycle 
duration sinusoidal phases of frequency 0.75 Hz 
traversing the VASAl1 array with slowness of 5.0, 4.8, 
4.6, 4.4, 4.2, and 4.0 sec/deg at an azimuth of 140 
degrees. In order of decreasing slowness the phases 
arrive at Eatonia at 5, 15, 25, 35, 45, and 55 ‘seconds. 
Each of figures 2.3 and 2.4 shows the input signal at 
Eatonia and the simple stack or 'beam' for a slowness 
of 4.6 sec/deg beneath the Covespagram. The random noise 
was synthesized by generating a 'random! pulse of the 


a t - e 
cos wt where A and w are random Gaussian 


form Re? 
variables and a> 0 is a fixed decay constant, at each 
digital point. The resultant noise was taken to be the 
sum of such random pulses. The Covespa patterns in 
figures 2.3 and 2.4 consist of contours of isolated 
"hills'; the contour value at the base of each hill 

is 0.5 and the contour interval is.0.1. From ‘figures 
2.3 and 2.4, we see that the Covespa process has 
successfully extracted the signals at the required 
Slownesses and times; the 'brightness' of the phases 


in the Covespagrams diminishes as the signal to noise 


ratio decreases. The superiority of Covespa as a 
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Covespagram for six synthetic events traversing 
the 1974 VASA. ‘Signal to noise aniplitude racic mere 
451.20 (he. verurcar trace immediately below the Coves- 
pagram is the synthetic seismogram at Eatonia. The 
coherent pulses at 5.0, 15207 2520 79 eo oo Uo 
55.0 sec cross the array at slowness values of 5.0, 
4.8, 4.6, 4.4, 4.2, and 4.0 sec/deg respectively as 
can be seen from the Covespagram. Below the Eatonia 
seismogram the 'beam' or resultant trace from all six 


Stations stacked at p = 4.6 sec/deg is shown. 
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Figure 24 


Covespagram for six synthetic events traversing 
the 1974 VASA. Signal to noise amplitude ratio is l:l. 
The vertical trace below the Covespagram is the 
synthetic seismogram at Eatonia. The coherent pulses 
at 5.0, 15.0, 25.0, 35.0, 45.0, and 55.0 seconds crocs 
the array at slowness values of 5.0, 4.8, 4.6, 4.4, 4.2, 
and 4.0 sec/deg respectively. Below the Eatonia seismo- 
gram the 'beam' or resultant trace from all six stations 


stacked at p = 4.6 sec/deg is shown. 
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resolving tool over simple stacking is clearly illustrated 
in these two diagrams; sole examination of the 'beam' in 
these figures would lead one to the incorrect conclusion 
that there is considerable energy arriving at a slowness 
of 4.6 sec/deg at times other than the true time of 25 
seconds. The Covespagrams in figures 2.3 and 2.4 have 
been calculated assuming the true azimuth of the incoming 
synthetic waves. Calculations assuming other azimuths 
reveal that coherency values, for a signal to noise ratio 
of 42d. for example, drop from above 0.9 to below 0.5 
for the coherent signals when the azimuth deviates from 
the true azimuth by 2 degrees. 

Figures 2.5, 2.6 and 2.7 show the Covespagrams 
for three events recorded by the 1974 VASA; the vertical 
traces of various stations have been used as input to 
the Covespa process. The Covespagram for a Salta Province, 
Argentina event (figure 2.5) exposes a P wave coda of 
over 60 seconds duration; similar P wave coda Covespa 
patterns have been observed by Gutowski (1974). The 
extended pattern can be interpreted as being the result 
of a finite time length source pulse convolved with 
source and receiver structure. Figure 2.6 shows the 
Covespagram from a Mongolian event as recorded at 
Eatonia, Maple Creek, and Suffield. Notice that the 
burst at about 9 seconds on the vertical traces is 
prominent on all stations. The onset of the event 


however is only clear on the Suffield record at about 
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Figure 2.5 

Covespagram for Salta Province, Argentina event 
(A = 83°, focal depth = 13 km, and magnitude = 5.5) 
as recorded at stations EAT, HAN, and SES. The vertical 


motion seismograms are shown below the Covespagram. 
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Figure 2.6 


Covespagram of event from Mongolia (A = 83°, 
focal depth = 33 km, and magnitude = 6.1) as recorded 
at EAT, MAP, and SES. The vertical motion seismograms 


are shown below the Covespagram. 
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FPigure.2.7. 


Covespagram of event from the Panama-Columbia 
border (A = 51°, focal depth = 5 km, and magnitude = 
5.4) as recorded at HAN, SES, and VUL. The vertical 


motion seismograms are shown below the Covespagram. 
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3 seconds. Nevertheless the Covespa technique does 
extract the onset as can be seen on the Covespagram. 
Again the P coda Covespa pattern is extended. The 
effect of a greater period pulse upon the Covespagram 
can be seen in figure 2.7 for an event from the Panama- 
Columbia border. The resulting Covespagram pattern 
becomes more bulbous as compared to the patterns of 

the previous examples. It is hoped that these Covespa 
examples together with those of Gutowski (1974) prove 
useful to future investigators employing the Covespa 


technique. 


The A(p), Slowness-Azimuth and t(p) Planes; Results from 


1974 VASA 


The 1974 version of VASA recorded teleseismic 
events (epicentral distance range of 30°< A < 100°) 
during wune ..ouly., and August, 1974, “Theyindividual: 
stations were operative during various overlapping periods 
within these three months. Forty-seven measurements of 
aT/dA and azimuth for the P phase were made using various 
elements of the VASA1 subarray. The various stations of 
the VASA2 subarray yielded sixty-eight measurements of 
adT/dA and azimuth for the P phase. Earth ellipticity 
and station elevation effects have been removed. Figure 
2.8 shows the epicenters of events for which VASA2 


observations of dT/dA and azimuth were made; almost all 
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Figure 2.8 


The earthquakes for which VASA2 observations 
of dT/dA and azimuth were made. The equidistant- 
azimuthal projection is centered at Hanna. The dashed 


circles are separated by 30° epicentral distance. 
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of the epicenters for VASAl1 observations are included 

on this figure. The results will be displayed here in 
the A(p), slowness-azimuth and tT(p) planes. [In all 
cases the observations will be compared to predictions 
based upon the assumption that the earth's lower mantle 
is spherically symmetrical with a velocity depth profile 


given by the Jeffreys-Bullen standard model. 


The A(p) plane 


The A(p) plane is of utmost importance in inves- 
tigations of the earth's lower mantle velocity structure. 
Depth ranges in the lower mantle for which velocity 
increases moderately and smoothly with depth are asso- 
Ciated with corresponding p ranges for which A increases 


smoothly with decrease in p. Rapid increases in velocity 


result in triplications in the A(p) plane and low velocity 


zones with sufficient magnitude are reflected by excur- 
sions of the A(p) curve for which A+ if the velocity 
reversal is smooth and for which A exhibits a finite 
jump increase if the velocity reversal is sharp. Now 
the Jeffreys-Bullen velocity model is such that velocity 
increases moderately and smoothly throughout the entire 
lower mantle; it is consistent with a chemically homo- 
geneous lower mantle. Departures in the lower mantle 
velocity profile from that given by Jeffreys and Bullen 


will have associated A(p) departures; the relationship 
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between A(p) and V(z) (Gerver and Markushevitch (1967) ) 
will then permit a refinement of the velocity-depth 
profile. This refinement in turn results in increased 
resolution of the lower mantle density, pressure, and 
temperature profiles via equations of state such as the 
Adams-Williamson relation. These considerations have 
been of prime importance in many studies. The following 
studies have revealed anomalous conditions in the lower 
Mantle; in all cases the 'anomalies' are manifested by 
the presence of large velocity gradients or discontinui- 
ties at various depths within the lower mantle. The 
large gradients are associated with ‘offsets' of the 
A(p) curve (rapid decreases in p over small A ranges) 

at various epicentral distances. 

Gutenberg (1958) suggested regions of depth 
900-1000 km and 1400-1500 km as being anomalous from 
amplitude studies. On measuring P and SH amplitude 
ratios, Vvedenskaya and Balakina (1959) concluded that 
there were anomalies at A = 38-42°, 51-53°, and 70°. 
Bugayevskii (1964) pointed to discontinuities in the 
travel-time distance curve at A = 36°=-37°, 51°-53°, and 
70°-73°. Carder (1964) reported that the travel-time 
distance curve could be represented as eight near 
straight line segments for the distance range 3°< A<102°. 
He states that the resultant lower mantle step model, 


with some modification and smoothing, would represent 
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conditions in the lower mantle more closely than the 
smooth Bullen model. Similarly Kanamori (1967) fitted 
the travel-time distance with three straight lines 

for the distance range 5°< A < 55°, His velocity model 
has rapid velocity increases at depths of about 150 km 
and 375 km. Chinnery and Toksoz (1967) studied the 
A(p) curve in the distance range 27° <1 90° sUuSing 
data from Lasa and noticed offsets in the A(p) curve 

at A = 35°, 53°, and 70° which could be correlated with 
velocity discontinuities at depths of about 800, 1300, 
and 2000 km. lLasa recordings of A(p) revealed anoma- 
lously high velocity gradients at depths of about 700, 
1200, and 1900 km (A = 35°, 52°, and 70°) an the study 
of Toksoz et al (1967); Rayleigh and Love wave disper- 
sion curves yielded ee wave velocity profiles with 
discontinuities at 350 and 700 km. Hales et al (1968), 
on the other hand, upon analyzing travel times to North 
American stations for 20° < A <.96° conclude that there 
is a discontinuity in the travel time curve at A = 24° 
and that there are no other major discontinuities up to 
A = 96°. Gopalakrishan (1969) used A(p) data for 

25° < A < 95° from the Gauribidancer (United Kingdom 
type) Seismic Array and compared his results with the 
Jeffreys-Bullen model; he found that his velocity model 
varied more slowly around ag fou km (fh = 30°). and more 


rapidly around Z = 1200 km (A =49°). Using measured p 
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values from 400 events recorded at Lasa, Chinnery (1969) 
found regions of anomalous velocity change near 700, 
1150, and 2000 km depth associated with anomalous por- 
tions in the A(p) for A = 32°-37°, A = 46°-48°, ana 

A = 65°-75°. On the basis of spectral amplitudes and 
travel times Archambeau et al (1969) found that P phases 
arriving at Western continental U.S. stations were 
affected by high velocity gradients at depths of 150, 
400, 650 and possibly 1000 km. The analysis of Johnson 
(1969), who used the Tonto Forest array, revealed 

OMe OMS eee )t. ete ome eee Ione ci OR Opt Ove tp) NATED 
81.5° corresponding to increased velocity gradients near 
Gepeicwotaco0,; 1 O00 aw! 230 to40 7) LOLO ande237 0) kim. 
Wright (1970) using the Warramunga array found regions 
of anomalously high gradients and also anomalously low 


gradient regions. Data from the four United Kingdom 


type arrays were combined by Corbishley (1970) who pointed 


to anomalous dT/dA features at 35-36°, 48-49°, 60°, 
68-70°, and 84-85° corresponding to possible variations 
in P-wave velocity gradient near depths of 850-900, 1200, 
1550, 1880-1900, and 2500 km respectively. Wu and Allen 
(1972), using stations of Wesson Observatory, noticed 
Oteront invites in) che: Ap) curve sai 52 ob 22. yands1G< 
butithey stated that it. would. be ditricult, to .conclude 
that these features were related to anomalous velocity 


gradients at the deepest portions of the ray paths. 


42 


170) oe i 


(way) Mer: 


00% wea oe 


tow endigeda, ,e aid te poate 

ae? ad casden se! ee wi vole 
aoeemiol Yo hen ann HAT apse aie 
as teaver Sayasve dheaade | 
Heo 4" 0 Gy (a5 62 ae A | ‘se Ue re 
wert Bina bere vaIDoisy eller dye jah 

7 Wy oe 

m1 OVES 06 Free weet ess boot — 
2Ibipet Dabo® YSTRE Arion? eS ops eurrou, Woe 
it Vietolenions Geib ban aE Cae pie pe 


; ee het 
soborksR hbethyt-avet e427 dearth sash pup lee ae 


eV 


bain com Gig (OVEL yvaidated: 800° penidaes: xeW ie 
,*0d , "Eb-oh OE BE Ye estuteeat bB\ Eb 2 om 


sPe 1 

enotssiszsy sldraeeqg oF rr oeeree "e+ be. & | 

A i ie 

OOS, ,000+088 to esq $60 tno hex tivokew 2 


AOD SA Bet ow yievidoeqner mand ae bns re 


he 


bao st ton »XTOd 20 ee cai %» noise | 


oy Bae 7S) ,°S2 $6 ave “wa mile: td on 


Pe aA he ee ‘ne 


sbrinntoo o3 MPD ERAS id hls, 28, an elle eT] | ea ay 
yoy ? 
ie vat ue mas 


ou 


viinoley avelemoas oo hotslox.s : 


Eee y 


* eriteq vex ‘eis to Be pnots10q | 


74 


Measurements of dT/dA from the Warramunga array 
(Wright and Cleary (1972)) revealed fairly low 
gradients at depths of 800 to 850 km, 1070 to 1110 kn, 
T1260 "tor 1330 Ki, © 7/50 “Pes kin, "and “2460 *Co"2600 km; 
high velocity gradients were found at depths of 1160 
to 1220 km, 2180 to 2370 km, and possibly 2700-2750 km. 
Vinnik et al (1972) studied dT/dA measurements at 
several arrays in the U.S.S.R. and found layers of high 
P wave velocity gradient near depths of 900, 1300, 1700, 
2000, and 2500 km. Hales and Herrin (1972) summarized 
the results “or Chrnnery “(1969 > "*Kattamor? (1967), Johnson 
(1969), and Corbishley (1970) and emphasized the features 
which these studies have in common. Wiggins et al (1973) 
used the A(p) data of Johnson (1967, 1969) and Wolfe 
(1969) in conjunction with an extremal type Herglotz- 
Wiechert inversion; their resulting velocity envelope 
reveals no fine structure in the mantle below about 
600 km depth. Kulh&nek and Brown (1974) inferred from 
UPSAS array data that the velocity distribution shows 
anomalous variations at depths of 750, 1500, 1800, 2300, 
anda=Z50@ km? | 
Thus there is much evidence that the P wave 
velocity profile of the lower mantle does deviate 
slightly from the smooth Jeffreys-Bullen model. Asso- 
Ciated with these small velocity excursions are anomalous 
perturbations of the A(p) curve. Also these perturba- 


tions are of a global nature since they have been reported 
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by various authors using different arrays on the surface 
of the earth. Figure 2.9 shows the A(p) results from 
the 1974 VASA. With event locations given by the PDE 
listings of the USGGS, all observations have been 
adjusted to a surface focus using the P wave velocity 
model of Haddon and Bullen (1969); it is felt that their 
“model is ideal for the purpose of focal depth corrections 
Since it is based upon free oscillation data and as such 
represents a global average. The distance coordinate of 
each point has been taken to be the distance from the 
equivalent surface focus to the center of the combina- 
tion of stations used for the calculation. The empirical 
solid curve is based upon numerical differentiation of 
the Jeffreys-Bullen travel-time table for a surface 
focus; a cubic-spline fit was made to a five point 
average smoothed version of the J-B tables using Fortran 
subroutine 'ICSICE' which is available in the Interna- 
tional Mathematical and Statistical Libraries. In 
figure 2.9 the small and large symbols correspond to 
VASA2 and VASA1 subarray observations respectively and 
the various symbols denote event azimuthal groups as 
indicated by the legend. Note that simultaneous calcu- 
lations of slowness and azimuth were made and that the 
observed slowness is plotted in figure 2.9 even if the 
measured azimuth may differ from the expected value; 


nevertheless such a A(p) plot should reveal 'radial' 
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Figure Poe. 


The 1974 VASA A(p) observations. The different 
symbols define 4 azimuthal groups; larger and smaller 
symbols refer to VASA1 and VASA2 observations respec- 
tively. All results have been reduced to surface 
focus; the solid curve is empirical and based upon 
numerical differentiation of the Jeffreys-Bullen travete 


time tables. 
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velocity anomalies if they exist. Notice that there is 
no particular A(p) signature characteristic of events 
from any azimuthal group or of determinations from 
either subarray. 

Now it is interesting to compare the 1974 VASA 
observations with the J-B curve and also with results 
reported by the various authors mentioned above, In 
particular let us focus our attention upon the major 
A(p) studies (those for which there were many observa- 
tions over the entire teleseismic range of 30°< A < 100°). 
These are the studies of Johnson (1969), Chinnery (1970), 
Wright (1970), and Corbishley (1970), which collectively 
cover a global distribution of observations in view of 
the various arrays used. These studies have exposed 
anomalously high velocity gradients at various depths 
and associated offsets in the A(p) at various epicentral 
distances; the results are summarized in Table 2.2. The 
most prominent offsets (those reported by the most 
authors) occur at distances of about A ~49°, A = 60°, 
and A = 70° and correspond to high velocity gradients 
near depths of about 1200 km, 1550 km, and 1900 km. 

Now the important question is - do the 1974 VASA A(p) 
observations exhibit offsets at the distances agreed 
upon by most of the authors cited in Table 2.2 (namely 
at 4249°, 60%, and /0?) an@wmi notpjiawhy? | Farstly, 


figure 2.9 reveals anomalously low p values (as compared 
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to the J-B curve) between about 46° and 55° and a sharp 
decrease in p values near about 50°. These features 

are similar to those associated with the 49° offset 
described by Johnson (1969), Chinnery (1970), Wright 
(1970), and Corbishley (1970). Hence, anomalously high 
velocity gradients near about 1200 km depth have possibly 
been detected by the 1974 VASA. Secondly, no statement 
can be made regarding possible anomalous conditions near 
1550 km depth in view of the absence of data points 
between 60° and 68°. Thirdly, it is inviting to associate 
the low p values in the distance range of about 72°< A< 76° 
with a p offset near 70° caused by anomalously high 
velocity gradients near 1900 km depth but the low density 
of observations in this distance range renders conclusions 
based upon the correlation rather weak. 

The 1974 VASA A(p) points describe a region of 
anomalously low p values in the distance range of about 
85° < A < 99° where obServed slownesses are as low as 
4.1 sec/deg. The J-B curve tapers off to a value of 
about 4.4 sec/deg beyond 90° and observations of the 
four authors mentioned above are in accordance with the 
J-B model. Thus if the anomalously low p values observed 
here in fact reflect anomalous velocity changes near the 
maximum depth of penetration of the associated rays 


(from a depth of about 2600 km to the core-mantle 


oe ae een 
ditpian ‘ (fen 


yidierog pase dash nd oost te 
tnenetste ort \ylbnoos2 ARAN b 
is9n enoitibnon svolemons. atyhoa in 
23ntog sisb 20 conta ont: to = car i 


"at >AR*SY Juods to spnset sonmetetb 
Heid ylevolsmons vd beaten, 0fse00 gontio. a 
ytierob wel oft tud dtqeb mt 000 ggen ewnbtbexe 
ecotaulonoo sxebass spars sonsdeth eind at | 
sow rsiie1 sOLseler zoo | 

to soiper 5 odizoaeb esniog iq)’ ABAY nae 4 | of bie | 
tuods to-spast eonaseth edz ai woulev g wot euc one pi 

es wol #8 9%6 agpesnwols pevtdedo e1reiiw is 4 ie = 

to suisy s oF te ategqss evita €-= ed? ab\ne 
eft to enolsevread bas *0C bioyed pab\soe bb 
ads dtiw soasbrocps nk ents cence hormkanem 4 ton 


bevxsedo eoulev q wol yYlavtoismons oy ea basal 


81 


boundary) then the anomalous velocities are certainly 
not a global phenomenon. Slowness values below 4.4 
sec/deg have also been observed by Kanasewich et al 
(1973) for rays which have their turning points in the 
vicinity of the core under Hawaii,by Hales et al (1968), 
and also by Davies and Sheppard (1972). The low p 
values observed here are associated with rays that 
bottom near the core under Hawaii and also near the 
core under other locations in the Pacific Ocean (these 
locations will be shown in the next section). The low 

p values for A> 85° could be produced by anomalously 
high velocities near the core mantle boundary; the 
increase in velocity needed to account for the A(p) 
observations depends upon the lateral extent of the 
anomaly (Green (1975)). It should be pointed out how- 
ever that both Kanasewich and Gutowski (1975) and Green 
(1975) report that low p values for A > 85° associated 
with near core rays arriving at such arrays as VASA and 
LASA from the west could be produced by anomalous Rocky 
Mountain upper mantle velocity distributions. Note 
also that Johnson (1969) has pointed out that the 
effect of the core should be considered when measure- 
ments of dadT/dA are made for the P phase for A >90°. 
This core effect only becomes significant when the 
difference between arrival times of the P and se phases 


is small compared to the period of the P phase energy. 
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The observed period of the 1974 VASA events for the 
distance range in question is about 1 sec and the J-B 
travel-time tables predict BrP oe travel-time differences 
greater ithan (0. 8 ‘sec for distances’ as great as’ ‘94°: 
hence few points will be effected by the core. Also 
Since the core correction when applied to the data of 
Johnson (1969) never exceeded about 0.08 sec/deg, it is 
felt that such a core correction would not significantly 
alter the 1974 VASA A(p) observations. 

A final question regarding the four major A(p) 
studies cited and the A(p) observations presented here 
remains. Although there is some agreement as to the 
location and nature of anomalies in the A(p) plane, why 
is there not total agreement? The answer to this mace 
tion is multifold. The entire teleseismic distance range 
is not well sampled - the density of observations as a 
function of epicentral distance varies from study to study. 
The degree of scatter of A(p) points about some average 
line is not the same for all studies and it is expected 
that purely radial velocity anomalies would be defined better 
by A(p) data sets with small scatter. The final and 
most important reason for the disparity is that a 
spherically symmetric velocity distribution is not an 
adequate description of actual conditions within the 
earth's mantle. Thus the results from each study are 


applicable to distinct regions within the earth defined 
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by the paths of rays between event epicenters and the 
array used and are not necessarily indicative of global 
conditions. Furthermore rays may suffer lateral dis- 
placements due to anomalous lateral velocity gradients 

in the course of their journey from source to receiver. 

In accordance later array studies of the 1970's 

included measurements of azimuth as well as slowness. 

The following section highlights the major results of 
these studies and includes a presentation and interpre- 
tation of the 1974 VASA observations in the slowness- 
Szamucho plane. ), Staal sit is hoped that the radial velocity 
anomalies revealed by the A(p) analyses mentioned are not 
viewed with extreme pessimism since they have been reported 
by many authors who collectively have studied a large 


portion of the earth. 


The Slowness-Azimuth Plane 


Figures 2.10 and 2.11 show the slowness-azimuth 
results for the 1974 VASAl1 and VASA2 subarrays respec- 
tively. In these figures there is one arrow for each 
seismic event with the tail of each arrow representing 
the slowness and azimuth predicted by the J-B tables 
for a ray originating at the uscGs location and arriv- 
ing at Suffield (for the VASA1 observations) and Hanna 
(for the VASA2 observations). The head of each arrow 


represents the observed slowness and azimuth. The 
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Figure 2.10 


Slowness-azimuth ‘array diagram' for the 1974 
VASA1 observations. The solid circles represent 
constant slowness values as indicated. The inner 
dashed circle represents the lowest possible J-B 
Slowness for the P phase (4.4 sec/deg). Each arrow 
represents one seismic event with the tail and head 
being at the expected and observed locations respec- 
tively. The projection is centered at the Suffield 


station. 
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Figure 2.11 


Slowness-azimuth ‘array diagram' for the 1974 VASA 
observations. The solid circles represent constant slow- 
ness values as indicated. The inner dashed circle 
represents the lowest possible J-B slowness for the P 
phase (4.4 sec/deg). Each arrow represents one seismic 
évent with the tail and head being at the expected and 
observed locations respectively. The projection is 


centered at the Hanna station. 
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length of the arrows are quite substantial (the average 
length for the VASAl1 subarray is about 0.31 sec/deg) 

and hence the immediate problem which arises is the 
identification of sources of the mislocations. The 

path of a P ray may be perturbed by anomalous velocity 
conditions near the source, the deepest portions of the 
ray path, and near the crustal and upper mantle region 
beneath the array. Generally it is accepted that the 
source region is the least likely cause of slowness and 
azimuth perturbations since the cone of rays which crosses 
an array, when traced back to a teleseismic source, is 
very narrow. Note however that source effects cannot be 
entirely ruled out (Davies and Sheppard (1972), for 
example). Past investigations of similar data have led 
to divided opinions as to the location of velocity anoma- 
lies which cause unexpected observed values of the vector 
ray parameter Pp (recall from Chapter 1 that a slowness 
and azimuth calculation is equivalent to a measurement 

of a 2 element vector ray parameter). Authors who have 
associated deep anomaly sources with teleseismic p 
deviations recorded at various arrays include Manchee 


and Weichert (1968) - Yellowknife array (YKA), Davies 


and Sheppard (1972) - LASA, Weichert (1972) - YKA, 
Kanasewich et al (1973) - VASA, Powell (1975) - Hanford 
array, and Wright and Lyons (1975) - YKA. On the other 


hand, vector ray parameter mislocations have been 
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attributed mostly to lateral inhomogeneities in the crust 
and upper mantle beneath arrays by Otsuka (1966a) - 
California array, Iyer (1971) - LASA, Naponen (1974) - 
Norsar, Hagfors, and Helsinki arrays, Wright et al (1974) - 
Warramunga array (WRA), Capon (1974) - LASA, Okal and 
Kuster (1975) - Tahiti and Rangiroa arrays, and Berteussen 
(1976) - NORSAR. Brown (1973) attributed p dislocations 
shown by the Uppsala array to both receiver and deep 
mantle effects. Also Engdahl and Felix (1971) and Aki 

et al (1976) have attributed anomalies in absolute teaven 
time of various body wave phases arriving at LASA from 
teleseismic events to lateral inhomogeneities in the 

crust and upper mantle under Montana. Julian and Sengupta 
(1973) and Sengupta and Toksoz (1976) studied body wave 
phase times on a world-wide basis and concluded that un- 
expected results could be attributed to lateral inhomo- 
geneities within the top 500 km and also within a region 
extending from a depth of about 2600 km to the core- 
mantle boundary. 

Thus from the above results it can be seen that 
caution should be exercised in the interpretation of 
‘array’ diagrams such as those shown in figures 2.10 and 
2.11. There are features other than the effects of 
velocity structure which give rise to unexpected values 
of p. Otsuka (1966) shows that the assumption of a 


plane wavefront with constant velocity introduces maximum 
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"errors' in time of about 0.2 sec for the case of actual 
teleseismic P phases impingent upon a 2° aperture array. 
Similarly it can be shown that the above assumptions 
introduce an error 6p such that || sp] | < 0.1 sec/deg 
for this experiment. Next it can be assumed that the 
errors in epicentral location given by the USCGS can 

be ignored (see Otsuka (1966) and Davies and Sheppard 
(1972) for example). Similarly the assumption that the 


surtield and Hanna Sites identify the 'center"® of the 


array for VASAl and VASA2 observations respectively intro- 


duces a negligible error. Also, as has been mentioned, 
earth ellipticity and elevation corrections have been 
applied. The sampling interval of the vertical motion 
is 80 ms and it has been pointed out that a cross 
correlation technigue provides time resolution of about 
20 ms. A generous estimate of the combined effects of 
the finite sampling rate and the Covespa process is about 
0.03 sec/deg for the combinations of stations used in 
this study. The value 0.1 sec/deg represents a maximum 
constant velocity-plane wave error and generally the 
error is substantially less than this; thus considering 
the intrinsic time resolution and calculation factors, 
it can be concluded that a realistic estimate of errors 
due to non-structural sources is || Sp] | Hi Ol Bac 


Now let Pr be the expected vector ray parameter, 


Pp 


Rowe the observed value, and define SPops by OD = 
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Pr 7 Pobs* As can be seen from figures 2.10 and 2.11 


the value ||ép 


aye is very often substantially greater 


than 0.1 sec/deg; thus anomalous velocity structures 
have affected the 1974 VASA observations. In an attempt 
to delineate the velocity deviations responsible for 


the large values of ||6p it is necessary to consi- 


Pops! | 
der the possibility of a 'deterministic' near receiver 
effect. The effect of dipping interfaces under an array 
upon slowness and azimuth calculations has been discussed 
by Niazi (1966), Manchee and Weichert (1968), and Zengeni 
(1970). Basically a displacement, 6Pei eb, which is 
approximately constant, is introduced; thus all observa- 
tions of Pp would contain a common component SBa, pee ifa 
series of dipping interfaces defined the crustal struc- 
ture under an array. In particular if there is only one 


dipping interface, such as the Moho say, the vector 


CPs ager is simply related to the strike and dip of the 


marker (Manchee and Weichert (1968)). Accordingly Capon 
(1974) has averaged the individual orthogonal components 
(North-South and East-West) of Pp observations recorded 


at LASA to determine the quantity 6p he then attri- 


=erust™ 


butes this constant component to near receiver determinis- 
tic structure and bases further interpretations of 


velocity anomalies upon a new p data set from which the 


bias has been removed. Such a calculation shows that 


SDs has magnitudes of 0.08 sec/deg and 0.15 sec/deg 
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for the VASA1 and VASA2 observations respectively. If 
a dipping Moho under VASA is responsible, then if the 
average crustal velocity is taken to be 6.5 km/sec and 
the Moho velocity is 8 km/sec the discontinuity has about 
1° and 2° of dip under the VASAl1 and VASA2 subarrays 
respectively. In both cases the direction of dip is to 
the north-east. It is felt however that the 1974 VASA Pp 
data set does not constitute a sufficient azimuthal 
sampling for such an analysis to be reliable and thus 
the observed 'bias' has not been removed. Results are 
given above to facilitate comparison with future studies. 
It is reasonable to describe the lengths and direc- 
tions of the arrows in figures 2.10 and 2.11 as being 
random in nature (casual inspection will verify this). 
For each event deviations of the travel time, 6t, at each 
station from the expected value give rise to the finite 
length arrows. The departures in travel time are in turn 
due to deviations of the velocity distribution from the 
standard spherically symmetrical J-B model. Now the 
slowness and azimuth of each event were determined using 
a particular combination of stations. Thus, in accordance 
with the results of Chapter 1, for each Pp calculation 
there will be associated expected deviations ép which 
depend upon the size of the travel-time perturbations, 6ét, 
and the combination of stations used. Figure 2.12 illus- 


trates this relationship; the abscissa is scaled such that 
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Figure® 2). LZ 


Diagram showing the ‘size' of the departures of 
1974 VASA slowness-azimuth measurements from quantities 
which would be predicted from the listed location and 
a J-B earth model. Each symbol represents one seismic 
event with the various symbols defining azimuthal groups 
and large and small symbols referring to VASA1 and 
VASA2 observations respectively (see figure 2,9). The 
observed deviation is given by the arrow lengths in 
figures 2.10 and 2.11. The expected deviation is the 
average ‘error' Opa (see Chapter 1) for random travel- 
time errors of 0.5 sec. Dislocations attributable to 
only non-structural sources have associated points which 
fall below line AB. LinesSAC and AD define expected 
deviations associated with random travel-time deviations 
due to structural effects of 0.2 and 0.5 sec respectively 


with allowance for non-structural effects. 
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it refers to the ‘average’ or expected error SPa (see 
Chapter 1) for random deviations in travel time of 0.5 
sec, the ordinate gives the lengths of the arrows in 


figures2.10 and 2.11 (that is the value ||6ép 


Penel ins and 


the points are plotted with the same convention as that 
for figure 2.9 (the various symbols refer to different 
azimuths of ray approach and large and small symbols 
refer to VASA1 and VASA2 observations respectively). 

The line AB is drawn (parallel to the expected deviation 
axis and intersecting the observed deviation axis at 


0.1 sec/deg) to represent the upper limit of contribu- 


tions due to non-structural sources. Thus, the deviation 


depicted by any point in figure 2.12 which falls below 
the line AB cannot be attributed to anomalous velocity 
conditions. The line AC is drawn such that departures 
depicted by points which fall within the region between 
lines AB and AC can be attributed to random travel-time 
deviations of about 0.2 sec arising from anomalous velo- 
city structure. Similarly the line AD is an upper 

bound for deviations due to structurally related travel- 
time departures of about 0.5 sec. The bounds given by 
lines AC and AD and the values of 0.2 sec and 0.5 sec 
associated with them are not intended to be 'sharp'; 

if an observation involving N stations has associated 
unexpected travel-time departures due to velocity struc- 


ture at the stations given by oti, tee pe ose Ne Glen ae 
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(Z8t2/n) * = 0.2 sec the ordinate of the line AC depicts 
the expected departure SP (see Chapter 1), which then 
depends upon the configuration of stations used,plus the 
constant maximum non-structural factor of 0.1 sec/deg. 
Similarly. if (26t2/n) 7 = 0.5 sec then the ordinate of the ; 
line AD depicts SPy plus the constant factor 0.1 sec/dec. 
Recall from Chapter 1 that observational deviations given 
by the form (26t4/N )% generally result in departures dp 
which describe ellipses. The average quantity Spy is 
considered here since the combinations of stations used 
were 'nearly' symmetrical and there was no one combina- 
tion employed a sufficient number of times such that an 
"elliptical' analysis would be much more appropriate. 
Notice also that figure 2.12 is based upon a least-squares 
type computation of slowness and azimuth whereas the 

1974 VASA observations result from Covespa computations. 
This computational disparity is not serious since a 
'Covespa plane wave' will closely approximate a best fit 
least squares plane wave. It is felt that the presenta- 
tion of figure 2.12 is a reliable and useful interpreta- 
tional aid. 

Now let us consider the effect that near receiver 
structure can have upon the results shown in figures 
2.10, 2.11, and 2.12. If we take some typical P wave 
crustal and upper mantle model, that of Massé (1973) say, 


then it is possible to determine the perturbation of 
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near vertical P wave travel times resulting from lateral 
changes in velocity beneath an array (the significance of 
the Massé (1972) model here is no more than that it is 
typical of models associated with stable continental 
regions. and if any other model were chosen then very 
Similar results would be obtained). Ten km undulations 
of the Conrad, Moho, top and bottom of the asthenopheric 
low velocity zone, and the 430 km discontinuities pro- 
duce travel-time deviations of about 0.20, 0.17, 0.02, 
0.07, and 0.10 sec respectively. Alternatively one could 
express lateral structural changes in terms of percent 
deviations. Thus, again using the Massé (1973) model 

as a basis, a 3% change in P wave velocity from the 
surface to a depth of 150 km produces travel-time changes 
of about 0.55 sec. If the magnitude of lateral velocity 
changes under VASA is ie described above then from 
figure 2.12 it can be seen that the 1974 VASA P obser- 
vations would be a poor indication of possible deep 
anomalies. It is unlikely however that near receiver 
structure is completely responsible for all of the 
observed anomalies. Firstly from figures2.10 and 2.11 

it can be seen that the pattern of anomalies given by 

the distribution of arrows for the VASA1 and VASA2 sub- 
arrays are strikingly similar - if near receiver structure 
was responsible for this pattern then the inhomogeneities 


underlying VASA1 and VASA2 would have to be nearly 
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identical and this is highly unlikely. Secondly there 
are regions in the slowness~-azimuth diagrams where the 
arrow directions and magnitudes change rapidly for small 
changes in epicentral location - this phenomenon is 
associated with velocity anomalies at deep portions of 
ray paths (Davies and Sheppard (1972)). Thirdly the 
extended Covespagram patterns shown in the previous 
section are typical of the 1974 VASA data and indicate 
that the velocity structure beneath individual stations 
is relatively uniform. 

A generous asSignment of travel-time deviations 
due to near receiver structure would be about 0.2 sec. 
In this case the points in figure 2.12 which lie above 
line AC are caused by deep anomalous conditions. 
Generally the slowness and azimuth of a given ray are 
most sensitive to anomalous velocity gradients, both 
lateral and vertical, near the deepest point of pene- 
tration. If we assume that observations for which values 


of ||ép fall above line AC reflect conditions near 


Pobs| | 
the associated ray bottoming points then several interest- 
ing statements can be made. The 'deep anomalies' then 

are associated with four major categories of SPops? 

(i) Low slowness observations near an epicentral 
distance of about 50°. Here SPops is largely radial in 


nature and points 'inward' on the array diagrams; the 


arrows associated with these anomalies shown in figures 
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2.10 and 2.11 have their tails between about 7.5 and 

7.8 sec/deg. This category has been discussed in the 
previous section where it was stated that the low 
slowness values could be associated with rapid increases 


in velocity gradient near about 1200 km depth. 


(ii) Low slawness values for some particular rays which 
bottom near the core-mantle boundary. In figure 2.10 
(VASA1) the arrows associated with these rays have their 
tails at slowness values less than 4.8 sec/deg at 
azimuths of about 177° and 310°. In figure 2.11 

(VASA2) the arrows have their tails at slowness values 
less than 4.8 sec/deg at azimuths of about 177°, 310°, 
and 240°. Again the anomalous arrows are mostly radial 
in nature and point inward. Notice that the group of 
VASA2 arrows at an azimuth of about 240° and between 
Sslownesses of about 4.4 to 5.0 sec/deg do not all 
conform; that is, some point radially inward while others 
depict near normal conditions. The bottoming points for 
the associated rays have surface projections near Hawaii 
and Kanasewich and Gutowski (1974) have noticed a similar 
complex anomalous pattern associated with this region. 

It has been pointed out that one cannot conclude that 


these types of 6p observations prove the existence 


=obs 
of lateral inhomogeneities near the core-mantle boundary 
especially when they are based upon rays arriving at 


VASA from the west. Notice, however, that this phenomenon 
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is also observed for rays approaching VASA at an azimuth 


On about. L772, 


(iii) A complex pattern of rapidly changing vectors 


Spans with a dominance of dp 


Pia vectors showing large 


azimuthal anomalies is the state which describes a 
series of observations centered be an azimuth of about 
138° and between slowness values of about 4.9 to 6.2 
sec/deg. This pattern which is similar for VASA1 and 
VASA2 can be seen in figures 2.10 and 2.11. The rays 
associated with this group bottom in a region between 
depths of about 1900 to 2600 km which has a surface 
projection near the Caribbean. It is interesting that 
in a study of seismic body wave phases recorded by the 
Yellowknife array, Wright and Lyons (1975) noticed 
rapidly changing azimuthal anomalies (from -4.3° to 
+1.4°) associated with rays bottoming below a depth 

of 2600 km under the Caribbean - they postulated the 
existence of strong lateral velocity gradients in the 
deep mantle under the Caribbean. Also, on the basis of 
travel-time residuals, Jordan and Lynn (1974) showed 
that anomalously high velocities exist in at least part 
of the depth range between 600 and 1400 km (and possibly 
deeper) under the Caribbean. The regions in question 
sampled by Wright and Lyons (1975), Jordan and Lynn 


(1974), and this study are all mutually disjoint but 


the results do suggest that the entire lower mantle below 
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the Caribbean and adjacent areas is in an anomalous 


state. 
(iv) Other spurious SPobs observations. This category 
includes the observations for which l}Spope! | falls 


above line AC in figure 2.12 other than those which 

fall in the above categories (i), (ii), and (iii). In 
general these observations constitute a low density 
sampling for any one region and they are scattered in 
the array diagrams. It is difficult to make comparisons 
of these 'anomalies' with those revealed by other 


studies because of their spurious nature. 


Figure 2.13 is an azimuthal great circle projection cen- 
tered at the Hanna station; the symbols in this figure 
depict the surface projection of the points of deepest 
penetration of rays for which VASA2 determinations of Pp 
were made. The general descriptions of SPobs refer to 
rays associated with bottoming points enclosed by the 
nearby boundaries and are based upon the observations 

of categories (i), (ii) and (iii) above. VASA2 mid- 
points are shown but the comments apply to the total 
array. Note that it is not necessarily intended that 
the anomalous observations are due to velocity structure 
near these points since, as has been mentioned, if near 


receiver lateral inhomogeneities are severe then no 


such claim can be made. Nevertheless, in view of the 
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Hanna showing the surface projections of the bottom- 
ing points of rays for which VASA2 slowness and 
azimuth measurements were made. The comments describe 
the nature of slowness-azimuth deviations associated 
with these rays and are based upon the VASAl and VASA2 


observations taken collectively. 
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fact that the interpretation of slowness and azimuth 
anomalies in terms of deep velocity structure is in 
accordance with results from many other studies it is 
felt that it should be considered as a serious possi- 
Diljty.. The VASA array is’ still inithe infant “stages 
compared to other arrays (LASA for example) and it is 
hoped that possible future VASA studies will provide 
detailed information regarding the near receiver struc- 
ture thereby rendering the results presented here more 
meaningful. Techniques described by Aki et al (1976) 
and Capon (1974) and the inherent mobility of VASA are 


encouraging indications that this goal can be attained. 


The tip Plane 


It has been stated that an exact interpretation of 
the p observations shown in figures 2.10 and 2.11 is 
not possible at this stage. Still it is worthwhile to 
consider the 1974 VASA observations in yet another plane - 
the t(p) vs p plane. Recently the importance of the 
function T(p) = Eee -p A{p) , where T is travel time, 
and its relationship to spherically symmetrical velocity 
depth profiles have been demonstrated by Bessonova et al 
(1974), Kennet (1976), and Bessonova et al (1976). By 


introducing the parameter of absolute time, T(p), into 


the previously discussed 1974 VASA p observations a new 
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Se6t of data’ t (p 


ae ) may be acquired where Pobs 1S 


obs 
the observed slowness and Tt (p ) is the observed 
bs ‘obs 


O 
value of tT. This data set is representative of the 
earth's ‘average’ spherically symmetric properties and 
as was the case for the discussion of A(p), azimuthal 
information has been discarded. Figures 2.14 a and b 
show the t(p) observations for the 1974 VASA; the p 
axis has been normalized such that a value p = 1 
corresponds to a ‘horizontal’ ray and an earth surface 
velocity of 6 km/sec. The actual slowness in sec/deg 


may be recovered by the relation p = (px6371x27) / 


actual 
(360x6). As in past diagrams large and small symbols 
refer to VASA1 and VASA2 calculations respectively and 
the various symbols identify event azimuth (see figure 
2.9). T(p) has been taken to be the average of the 
observed travel times for the combination of stations 
used corrected to surface focus and A(p) has been taken 
to be the distance between the event epicenter and the 
center of the combination of stations used, again 
corrected to surface focus. As was the case for the 

A(p) study the Haddon-Bullen (1969) P wave velocity model 
was used for surface focus corrections. The empirical 
line in figure 2.14 is based upon numerical differentia- 
tion of the J-B P wave tables for a surface focus. The 
technique of differentiation has been described in the 


section concerning the A(p) plane. Notice that the 
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Figure 2.14 a and b 


The 1974 VASA t(p) observations. The various 
symbols refer to azimuthal groups with larger and 
smaller symbols associated with VASAl and VASA2 obser- 
vations respectively (see figure 2.9), The empirical 
solid curve is based upon numerical differentiation of 
the Jeffreys - Bullen travel-time tables of P waves for 
surface focus. The slowness axis, p, has been nor- 
malized assuming a surface velocity of 6 km/sec. Actual 
slowness values in sec/deg may be recovered using the 


relation p (p x 6371 x 27) / (360 x 6). 


actual 
The set of axes applicable to the data shown are 


in accordance with the inequality dt/dp < 0. 
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calculated tT(p) values are decidedly lower than would 
be predicted by the J-B tables especially in the p 
region of 0.30 to 0.48 (5.50 to 8.9 sec/deg). Also the 
J-B curve proceeds to no lower p values than about 

0.24 (or 4.4 sec/deg) whereas the observed values extend 
to p values less than this. This anomalous feature has 
been discussed in the previous sections. Now Hales et 
al (1968) have provided an expression between T and A 
based on teleseismic observations. The relation which 
assumes a source upper mantle velocity structure asso- 
ciated with tectonically active regions and a receiver 
upper mantle structure associated with stable continen- 


tal regions is, 


(A) = 72.774 10.9D10A - 3.2087 x 10°07 A~ = 2.003% 1077 A> 


G2ii2)) 


where A> 30° and T is given in seconds. In view of the 
location of events (figure 2.8) and the array used in 
this study it is felt that the above expression should 
be comparable to the 1974 VASA observations. Using the 
fact that le = -A(p) and hence that t(p) is a mono- 
tonically decreasing function (see Chapter 3) one can 
construct an empirical t(p) curve based upon the T(A) 
relation above given by Hales et al (1968). Figures 
2.15 a and b show the same tT(p) observations as those 


in figures 2.14 a and b but the empirical curve is now 


109 


. Bbnedxe aenutsv purnesad atts asoradw ( 
eet suvtset evolsmons edt ati a 
te a9isH wot .anoliose ail els ak be 
4 bos T nsewsted cotseoxxe ae babivona evad ( 
fioidw noitsf[et ofT .anotdsvxsedo ‘ous | 


~oaas owwjouise ytioolev saan voto pean’) 


xevieosi 5 Sis” enolpex svidos rhs: 


-fentisznoo sidate adiw ‘Baise a toux: | 
# vat enoipe: | 


€,°-or w £00.58 Rae € = AOLSe, OL + 8% “ST = (a) iy . 
we 


(S.&) | : r 


eit to woiv nl .ebaesse ai aavip at . 5s. °0€ <A ; 
ni bees yerrs sd¢ bas (8.8 smupi2) eineve Fo aoitsooL : ; 
binole nolzeetqxe evods ont tad¢ ist ef 3 ih aids 
att opaieu .enoissyvreado ASAY BYGI saz oF olden : 
~onom & et (q) 7 tedd Sono Bas (q}4- = juts $088 
mso ene (€ tejqsiD sea) soliton? sities eiteokned 
(A)T edd noqu beesd eviso ape fsoixiqme as Er rade 
Titi bde .(83@L) Is do stow * ee ae | | 


V4 i 


) a 4 BT 
" iy ha ay 
Ta oa Ty ie ee | 
: Perks te $i a 
ad Wo), (9 in 

7 j ts i 


oan a ie me hel Ale inj nt 


ae 
mg Oe Mame Aer ee 


o. % | 0.40 
# | 


Figure 2.15 a and b 


The 1974 VASA t(p) observations. All descriptions 
are the same as those for figures 2.14 a and b except 
that the solid empirical curve here is based upon the 


T(A) relation (equation 2.2)) of Hales et al (1968). 
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based upon the above T(A) relation; there is a striking 
Similarity between the observations and this empirical 
curve. Unlike the case farthe J-B comparison, points 

no longer fall decidedly away from the solid line and 
the Hales et al (1968) t(p) curve extends to a p value 
ee as about 0.227 (4.20 sec/deg) approaching the low 
values of the observations. The Hales et al (1968) 
velocity depth profile is consistent with their T(A) 
relation to within 0.01 sec. The velocities are always 
greater than the J-B model with a maximum deviation of 
about 0.1 km/sec occurring at a depth of about 1200 km 
and substantially smaller deviations throughout the 

rest of the lower mantle. Thus the 1974 VASA t(p) ob- 
servations are in accordance with an 'overall' spheri- 
cally symmetric P wave velocity profile which is slightly 
greater than the fae of Jeffreys and Bullen. The 
empirical T(A) relation of Hales et al (1968) should 

be considered a serious candidate for a basis of com- 


parison for future array P wave studies. 


Conclusions 


The Covespa technique proved to be useful for the 
slowness and azimuth calculations of teleseismic P waves 
crossing the 1974 Variable Aperture Seismic Array. 
Deviations of the slowness and azimuth observations from 


quantities which are in accordance with a spherically 
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symmetrical earth for which the kinematic properties of 
teleseismic P waves are defined by the Jeffreys-Bullen 
Seismological Tables are in accordance with the existence 
of lateral P wave velocity inhomogeneities along the 
paths of rays between event locations and VASA. If 
lateral velocity changes in the crust and upper mantle 
under VASA are severe then the slowness and azimuth 
observations would be a poor indication of possible deep 
mantle anomalies. However preliminary indications are 
that this is not the case. If vector ray parameter 
deviations are associated with anomalous velocity con- 
ditions at the ray bottoming points (the region to which 
vector ray parameters are most sensitive) then low slow- 
ness values near A y 50° are in accordance with anoma- 
lously high velocity-depth gradients at a depth of about 
1200 km, low slowness values for A > 85° are consistent 
with velocity inhomogeneities near the core-mantle 
boundary under some areas of the Pacific Ocean, and 
azimuthal deviations associated with rays passing under 
a region near the Caribbean point to anomalous velocity 
conditions predominated by lateral velocity gradients at 
depths between 1900 and 2600 km under that region. Also 
the 1974 VASA slowness observations can be expected for 
rays emerging from tectonically active regions, passing 
through an 'average' spherically symmetric lower mantle 


for which the P wave velocity is slightly greater than 
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the Jeffreys-Bullen designation, and impinging upon a 


seismic array in a stable continental region. 
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CHAPTER 3 


INVERSION OF TRAVEL-TIME DATA USING THE TAU METHOD 


Introduction 


In past linear seismic arrays have been used widely 
for the investigation of crustal and upper mantle struc- 
ture. Travel-time distance curves have been inverted to 
obtain spherically rept hao velocity depth functions 
using various inversion technigues. Inversion techniques 
may be divided into two major categories; indirect methods 
(for example, Backus and Gilbert (1970)) and direct methods 
such as the Herglotz-Wiechert techniques. Attention here 
will be focused upon powerful technique of seismic travel- 


time inversion, the Tau method, which was developed by 


Bessonova et al (1974). The method makes use of the 
eb igtod ae key ¢) 
T(p)) =p) -- Sp x ep) (Sens) 


where p = daT/dX is the ray parameter, X is half the epi- 
central distance, and T is hakf the travel time from source 
to receiver. It was shown that it is possible to map 
lamitsiofitherfunction) tiip)eimto limits» forethe velocity 
depth function V(y). Thus with t(p) estimated from T(X) 
observations an envelope may be obtained in the V(y) plane 


which contains all possible velocity depth profiles that 
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are consistent with the data. Prior to the development 
of the Tau method ene Herglotz-Wiechert method of inver- 
Sion, which maps the function X(p) into V(y), was widely 
aSeaqhrorsdt@rect “inversion. “in its cYassical”form the 
Herglotz-Wiechert method was restricted to velocity- 
depth profiles for which there are no velocity reversals 
(Herglotz’ (1907) and Wiechert:(1910)). Gerver and 
Markushevich (1966) overcame this serious restriction 

by developing a 'Herglotz-Wiechert type' expression which 
accounted for the presence of low velocity zones. But 
this inversion scheme still suffered from a serious set- 
baék*sinee it implicitly assumed that the function’ X (p) 
was known exactly. Accordingly McMechan and Wiggins 
(1972) extended the Gerver-Markushevich formula to an 
extremal inversion scheme which mapped limits in the X(p) 
plane into limits in the V(y) plane. The disadvantage of 
this mechod- as) thak.wt as ditficult,.to: put bounds on. the 
Function (sp) since X(p)>'@ when fa ‘Smooth velocity rever- 
sal is encountered. The function t(p), on the other hand, 
is well behaved since it decreases monotonically and 
responds to a low velocity zone by exhibiting a finite 
discontinuity. It is expected that» the Tau method of 
seismic inversion will become a standard tool in the in- 
terpretation of crustal, upper mantle, and teleseismic 
data. 


In this chapter a new method for obtaining bounds 
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for the* function’ tT (p), ‘using data “from Tirnear ‘arrays? is 
presented. For each branch of the travel-time curve, T 
observations are fitted to a family of second order 
polynomials in X. The families of curves are then mapped 
into the t(p) plane. The Tau method will be illustrated 
by inverting T(X) data recorded by the University of 
Alberta on Project Early Rise. A comparison of the results 
of this study with those of numerous other authors will 

be, given. , Also, it ss stressed that uncertainty in velo- 
city depth profiles is due partly to errors in the quanti- 
ties, TI, Xs. andjpenand at ius) also. affected by thesfact that 
observations Of I, “4, sends p are incomplete; thus the 
resolving power, as a function of separation of observa- 
tion points, of the Tau method is examined by inverting 
tT(p) envelopes calculated from exact velocity-depth 


Pac cCLoOrs - 


Mapping “t(p) Limits into Velocity-—Depth Limits 


The following analysis highlights the essential 
features of the Tau method described by Bessonova et al 
(1974) and also includes some interesting points which 
became apparent during various calculations. 

It is convenient to discuss the inversion problem 
in terms of a flat earth model; the flat earth transfor- 
mation makes this possible. Suppose we have a spherical 
' 


velocity function v(r), where "r” is the radial coordinate 
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from the center of the earth, and associated body wave 
observations T(p), A(p) and p = dT/dA where T(p) is 
the half travel time and A(p) is the half epicentral 


distance. The normalized earth flattening transformation 
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u(y) = Ves ___. | (3.2) 
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The normalized flat earth model defined by (3.2) is 
characterized by horizontal distance coordinate X(p), 
general depth coordinate, 'y"', and slowness depth profile 
u(y) where u(0)=1. Rays in this flat earth satisfy 
Snell's law, sin a(y) =pu(y) where a is the angle between 
the y axis and the Girection of the ray. The problem now 
reduces to mapping limits of the function tT(p) = T(p) - pX(p) 
into upper and lower bounds for the depth, Y(p), at which 
the ray with parameter p bottoms. It is extremely impor- 
tant to realize that the mathematical expressions required 
for inversion, and indeed. the properties of T(p), Xp), 
t(p) and other functions related to seismic ray propaga- 
tion, are based upon the following assumptions regarding 


the slowness-depth profile u(y) given by Gerver and 
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Markushevich (1972); 

(a) The function u(y) is positive and scaled so that 
u(0) =1. (Note that this is ensured by the transforma- 
THON 03.) Ju 

(b) The function u(y) is everywhere twice continuously 
differentiable, with the exception of a finite number of 
points at which it or its derivatives either do not exist 
or are discontinuous. 

tc) In every finite segment, u(y) is bounded, but on 
the entire semiaxis y > 0 it is not necessarily bounded. 
(d) There exists a finite number of waveguides (low 
velocity zones). If we define the function s(y) by 

s(y) = sup{u(y°?), y° « [0,y]} then the waveguides are 


intervals of the y axis that have the following proper- 


ties; 
C2) s(y) is constant in each of them. 
(ii) Each of them contains points where 
Gy) <7e. (Vs. 


(iii) Outside these intervals, u(y) = s(y). 


The assumptions (a), (b), (c), and (d) above clearly do 

not impose any serious limitations when one is dealing 

with practical spherically symmetric earth models but as 
has been mentioned they are needed for the mathematical 
development of the inversion problem. Now the low velocity 


zones are ordered by the index "i" and are described by, 
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ea depth to the top of the yen low velocity zone. 
Y; = depth to the bottom of the neu low velocity zone. 
qo = slowness at the top of the eas low velocity zone. 


Notice that the functions Y(p) and u(y) are mutual 
inverses outside the low velocity zones; thus flat earth 
bounds for Y(p) may be transformed into bounds for the 
spherical earth velocity function v(r) outside the wave- 
puidessusing (3,2)... lt is important (to realize that it 
is implicitly assumed that the travel-time curve is 
composed of unconverted body wave phases and that if 

a ray is reflected then it is totally reflected. A 


precise mathematical statement of this assumption is, 
Vito) Gabino ty, ounly): el by efor “pre (0g0))-.. 


Now the solution of the inverse problem is (Gerver 


and Markushevich (1966)), 
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Y(p) = aaa | X (q) Iq? = p7] dq 
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Let us consider the average depth of bottoming 


¥(m, ,m,) for a package of rays m SP<M,, 


e -1 /2 
mel 
From’ (3..3), and (3.4), ¥(m,,m,) can be written as 


¥ (m, ,m,) = o(m,,m,)- y(m,,m,) . (3.5) 


The function ¢(m,,m.) is related to the Herglotz-Wiechert 
portion (first term) of the right hand side of (3.3) and 


it is given by 


ay -1 [72 ie <8 
o(m,,m5) = 27 (m5 - m,) | (q)m,q (g -m,) dg 
uh a - 
+ | rq) [mya7* (a? = mj) #~ myq™* (a? - mb) Jas} (3.6) 
™m 


If no low velocity zones are present or if for every low 


velocity zone q;<m then ¥(m,,m,) = >(m,,m,). The term 


ie 


v(m, -m,) is related to the low velocity zones and is given 


by 
y(m,,m,) = Oa ), Bena.) mM.) (337) 
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arcos (m,/q; ) LG m, $4; <M, 
a(q,7m,,m,) = | 

arcos (m,/q,) - arcos(m,/q;) Aisa m, Sq; s ake 

(3.9) 

The summation in (3.7) is taken over all low velocity 
zones for which q3? m+ It is possible to show that 
B, (m,,m,) > 0. Also, when m5< qj, then there exists a 
value jy, my Say < M5, such that 


a 


B, (m,,m,) = Lous pang c4Ok dy (3.10) 
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_where 
ee 
Q.= {luty)-af]/lag- 71}. 


The expression (3.10) for B. (m,-m,) Ee. Dabrereu bak ay 
useful for two reasons. In practise we do not know qy 
and O. exactly. From the T(X) data however we can find 


limits for qe and Oo. (Bessonova et al (1974)) such that 
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Using (3.10), the upper bound, B.(m,,m,), such that 


B, (m,,m5) bs B, (m,,m,) is given by: 


B, (m,,m,) = om eke anand Lie Oo. tan", ] (371.45) 


where 


a 
an 2 - 2 2 Dene 
Thus B. (m,,m,) is maximal when the low velocity zone is 
rectangular and takes on the maximum possible velocity. 


th low 


Furthermore, the maximum thickness, Bag for the 7 
velocity .zone can be found using, the induction method 
described by Bessonova et al (1974). The lower bound, 


B,(m,,m5), such that B, (mj ,m,) > B. (m,,m5) is then given 


Dy. 
2 Be a 1 
B. (m,,m,) = G,la;-y ] [1-Q,"tan Q.] 63. 1:2) 
where 
a D =D -2 gic % 
and 


2,72 -2 
{os/he + q.} 


Thus, B. (m,,m5) is minimal when the low velocity zone is 
rectangular and as thick as possible. The form (3.10) 
for B; (m,,m5) is instrumental in the analytical deter- 
mination of the upper and lower bounds, (3.11) and (3.12) 
respectively, for B, (mj) ,m,). It is also useful for com- 


putational purposes. This can be seen by examining the 
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expression for B, (m),m,) given by (3.8); the second term 
on the right hand side of the expression involves an 
integration, with respect to the variable 'p', which 
cannot be performed explicitly. The expression (3.10) 
for B, (m,-m,) is exact for some unknown value of y such 
that my ORY, <m.,- 


exactly, approximations for the right hand sides of (3.10), 


However, even though we do not know y 


(3.11), and (3.12) may be™made by taking y= (m,+m,)/2. 
Of course, we expect the approximation to be good if the 
interval (m, ,m,) TS (smauL. I have performed several 
calculations comparing the 'exact' value (3.8) with the 
approximate value (3.10), where y = (m,;+m,)/2, FOrSrees 
tangular low velocity zones. The integration with respect 
to, p>, 1n (3.28), was, done numerically.” in ald cases, it 
was found that results from (3.8) and (3.10) compare 
favourably. 

We have upper and lower limits, B. (m,/m,) and 
B, (m,,m,) respectively, for the functions B.(m,,m,). Now 
suppose that we have upper and lower bounds for the func- 


tion t(p) ian the-ainterval (d,1) such that 
7 (pypPend (pyv< Vs tp} for OWS Se Dee Ls 


Then for m, >d there exist two functions, £(m,,mo) and 


1 
g(m,,m,), which bound >, (my -m5): 


£ (m, ,m,) <>, (m,,m,) s g(m,,m.) (23) 


rae Lt ae 
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ieee 4 cif 
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wom .( cits pM) G& anoisonv? sis rot .yisvisosgees tga ik) 
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bak Lg a ee tt | ie 


ue Bie oe ehé al 


eae ( sahepdnna. lati ) e 
"i i 


2 
# 


m 
2 
-l =1 — hs, 2he2 
£(m,,m5) = 27 (m5-m, ) | t(q)m,q (g -m)) dg 
m 


- 1 
2 ie, wbbeowt ec i@. Dhawan 
+ | t(q) (m,q (q -m) -mjq ~(q -m.,) )aq (3 214) 
mS 


and 


m 
1 
oe 


2 
-1 -1 ~ ro ee ae” 
g(m,,m,) = 27 (m5-m, ) | tT(q)m,q (g -m,) dq 
+ | t(q) (mq *(q -m3) = mq” > (q*-m5) Jaap. (305) 
m 


Now since Y(p) iS a monotonic non-decreasing func- 


tion we have, 


Y (m) 


1A 


¥ (m, ,m5) for m 


IV 
= 


: CSeeeG)) 


1A 
3 


Yecm) 2 ¥ (mj ,m,) for m 


BernoetoL Ge ase L)e (ato We la) ee 3. Wd eee et bo); feand 


(3.16) then, it is possible to find upper and lower bounds 


for the Tumction Y(p). ip’ the anterval “de pr<1... We assume 
ct Rg (seas T(p) fonvd: = pi. and that there are "SIL" 


non overlapping low velocity zones defined by, 


(05105095 /45-¥;-4;) i = 1,NL 


For simplicity we assume that the first low velocity zone 
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does not commence at the surface and also that the last 
low velocity zone does not commence at the deepest 

point of the observations. Thus qi < 1 and Wnt 7 ds 

dn the, absence of low velocity zones, the upper and 
lower estimates of aL te Y (p) and Y(p), are given 


Dye 


Y (p) inf (g(m,,p) : Shas m, < p) (Sey) 


I 


¥ (p) sup(f(p,m,) : eres ie (3.138) 


Now, let r(p) and s(p) be defined by: 
f(p,r(p)) =sup(f£(p,m,) : p<m, & Yr (35°19) 
g(s(p),p) =inf(g(m,,p) :d < m, < Dp) (3.20) 


In general, we calculate upper and lower bounds for Y(p) 


for p values Ee j= 1,N) where ie 8 When low 


itd 
velocity zones are present, the upper and lower estimates 


are given by: 


Y = .}- ; oe 382 
Y@,) g(s(p,) p,) cy B, (s(p,) P5) (30°20) 
eS | 
= ))- B. oe ye 3.22 
¥(p;) £(p,,r(p,)) pel ye B, (p, r(p;)) ( ) 
Peso fe 3 


In using (3.21)7and (3:22), it. 1s impontantito: be aware 
of the following guidelines. First of all, calculations 
using (3.21) and (3.22) are done in order of decreasing p 


values; that is in order of increasing "j" values. In 
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(3.21) the summation is taken over all low velocity zones 
for which q, > P;- Alsoyw vf therinverval (p,,r(p;)) con- 
tains a low velocity zone, [that is, if there exists a 
q, such that q, <q, < q, and al veel r(p.)l, then 


(3.22) is replaced by ees = ei ). Furthermore, 


a 
notice that the values B; (s(p;) -p5) are a priori unknown 
for non zero Oo; Since they require, from (3.12), knowledge 
of undetermined values h, - Thus, the following induction 
process is performed. In general, hy is the supremum of 
the difference ¥(p;) - ¥(p5), calculated for P5 in the 
region gq, < Pas qn where ¥(p;) and ¥(p,) are given by 
(352) antea 3.22) with allowance for ‘the first (=1) tow 
velocity zones. Thus, ae is determined from (3.21) and 


(3.22) assuming no low velocity zones. Then, in turn, 


hoha,---/hy are determined. Once all of the maximum 
thicknesses for the low velocity zones are given, ¥(p5) 
and Lee, are given unambiguously by (3.21) and (3.22). 
NetiGerthat thevwtuncii ons £(m,,m,) and g(m,,m5), 
given by (3.14) and (3.15) respectively, play an important 
role in the inversion process whether low velocity zones 
are Dr esention mot... thn general, when 1. (p)e and t(D). are 
given as piecewise polynomials, of up to second degree in 
p, then (3.14) and (3.15) may be integrated explicitly to 
give £(m,,m,) aS a piecewise smooth function of m. for 
fixed m, and g(m,,m.) aS a piecewise smooth function of 


1 


my for fixed M. « The problem then is to determine 
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sup (f(m,,m,) : m)<m, s 1) for fixed my 


Gus m, < m. } for fixed m4. In general, this, is. done. by 


sampling £(m,,m,) for various m. and g(m,,m.) for various 


and inf (g(m,,m,) : 


mM, - There is at least one situation, however, in which 
sampling need only be done at 'preferred' points. For 


example,..suppose..that t{p) = t.w(asconstant).and t(p)=T 


U L 


(a constant, Ty 7 Ty) for the interval qs pie G5: Then, 


for fixed m, <q), the maximum value of f(m,,m,), for m. 
restricted to the interval qd, £ ™ <£ do, WilLoccur.at 
either M5 =, Or M4=q>. This follows from the fact that, 


in this case for q, < Mm, S$ do, £(m,,m,) is doubly smooth 


3°£(m,,m5) ‘im 

and Tsemtonie oe Oe Similanly form @ixed mM. > Aor the 

m 
2 


minimum value of g(m,,m for m, restricted to the inter- 


2)s 1 
val qy SM, S$ dor Wield .OCcCuUm vat either m,= 4, Or M,=q>- 
According. y),. for Gq, <™ 5G, g(m,,m.) is doubly smooth 


3°g (my ,m9) 


and Oe. 
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Resolving Power as a Function of Station Spacing 


It has been shown, (Herglotz, 1907) and (Wiechert, 
1910) that in the absence of low velocity zones if X(p) is 
known exactly then V(y) iS unique and can be calculated. 
The presence of a low velocity zone introduces uncer- 
tainty in the velocity-depth function below the low velo- 
city zone. Gerver and Markushevich (1966) have obtained 


a theoretical expression which permits one to compute the 
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uncertainty in the velocity-depth function when any 
number of low velocity zones are present assuming that 
the functions X(p) and T(p) are completely determined. 
In practise, however, seismic observations are taken at 
discrete eee Xa, along the surface of the earth. 
Thus, even in the absence of low velocity zones, there 
will be uncertainty in the velocity depth profile due to 
the discrete nature of the observations. Davies and 
Chapman (1975) have discussed the discrete inversion 
problem, restricted to a travel-time branch within which 
neither triplications due to rapid increases in velocity 
nor discontinuities due to low velocity zones take place, 
in terms of the classical Herglotz-Wiechert method. 

It is possible to study the effect that discrete 
T(X) data has upon the resolution of the velocity depth 
profile using the Tau method. It is assumed that exact 
observations of T, X and p are obtained from N small 
linear arrays of seismometers with array spacing being 
AX and the array lengths <<AX. In general, a given array 
may not provide a set of observations (T,X, 0Pys i= 1,N) 
since Xs may be located within a shadow zone and also a 
given array may provide more than one set TT; 7X3 9Ps since 
x. may be located within a region of triplication. For 
this example the function t(p) will be due to the dis- 
crete nature of the Tt; (p;) observations. Since Se = 7X (Pp), 


and hence t(p) is monotonic decreasing, it is possible to 
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find functions Ty (P) and tT, (p) from the observations 
Tt, (p;) such that «(p) must satisfy, Ty (p) < ig Vee eS Ty (P) - 


For each interval (D544 7P,) where P;>P the upper and 


sie pl 


lower Tau bounds are given by, 
ie ane eee 


Ty (p) = Tt, (p,) : 


Clearly a smaller station spacing AX will result ina 
small t envelope defined by Tr (Pp) and Ty, (p) since there 
will be more observations tT; (p,). A narrower t(p) envelope 
will in turn result in a narrower velocity envelope. An 
example of this is shown in figures 3.l1aandb. Figure 
3.1 shows a hypothetical P wave velocity depth function. 
Using a direct ray program it is possible to calculate 
values T; (p,), X.(p;) and hence Tt; (p;) generated by such 
a structure using the Bullen (1963) ray integrals. 
Figures 3.l(a) and 3.1(b) show upper and lower velocity 
bounds calculated from theoretical observations based on 
the velocity model in figure 3.1 for station spacings of 
100 km, and 25 km respectively. The inversion routine 
used for mapping t(p) limits into velocity-depth limits 
is after Bessonova (1974) and has been described in the 
preceding section. A similar analysis may be performed 
for any hypothetical velocity depth structure and thus 
aid in the determination of station spacing required for 


any given desired resolution. 
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Upper and lower velocity-depth bounds resulting 
from Tau inversion. In parts (a) and (b) the symbols 
represent the upper and lower bounds assuming exact 
knowledge of t(p) at intervals of 100 and 25 km 


respectively based upon the solid velocity-depth curve. 
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Estimation of the Function t(p) from Real Data 


Before inverting travel-time data using the Tau 
method it is necessary to find estimates for the function 
tT (p) from T(X) observations. Bessonova et al (1974) 
formulate the problem of t(p) calculation using the 


Clairaut equation 


T(p) = p go + al- a : (3223) 


The set of particular solutions of (3.23) is given by 


Ceol). 2f rom the fact that p= dT Bessonova et al (1974) 


Axa! 
Show that 1(p) jis. thevsindgular solution? of) .(3623)); that 
is t(p) is equal to the value 1t(X,p) for which ou ieiP) = 0. 
Thus for a generally forward branch within which the slope 
p= is present T(p) will be equal to the maximum value 
of t(X,,p) calculated for all X.. For a generally 
receding branch T(p) will be the minimum value of T(X;,/P)- 
An example of this is shown in figure 3.2 for the forward 
branch consisting of rays which bottom within the first 
layer of the hypothetical velocity structure shown in 
figure 3.1. The peaks are very clear; in addition, as 
expected for a forward branch, the X coordinates of the 
peaks decrease with increasing p values. Bessonova et al 
(1974) use this regular behaviour of the peaks as an 
interpretational aid. 

The method of t(p) calculation described above is 


excellent when errors in T(X) data are very small. On the 
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Figure 3.2 


The function T(x; p) "for the “first forward trave.— 
time branch arising from the hypothetical velocity curve 


of figure 3.1. The p values are such that Psay > Pi 
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other hand, when even modest errors are present or when 
the elastic waves traverse layers with moderate hetero- 
geneities it becomes increasingly difficult to extract 
(ph) eucormataon using “Clairaut's relation” ay btgures 3.5 
Bhone the travel times and distances of the first branch 
from a seismic line recorded during Project Early Rise by 
Creve ersity of Alberta. 

The graph of T(X. Pp) for p values Py? B57 P37 Pa? Ps 
arising from this generally forward branch is shown in 
figure 3.4. Because of inhomogeneities in the layered 
media and possible errors in time, the regular behaviour 
which is exhibited by hypothetical example is non seen. 


In particular, for any fixed p value, there are several 
ox 


travel-time points of figure 3.3 is typical of long range 


X positions at which is zero. The scatter in the 
refraction surveys, and is in part due to local variations 
in crustal structure. Thus it is necessary to smooth T(X) 
in some manner in order to derive spherically symmetrical 
velocity depth functions. Kennett (1976) suggests that 
this may be done either by fitting a lightly smoothed 
spline.to the function tT (Xp) for fixed values of p or 

by interpolating each branch of the travel-time curve with 
a smoothed spline which is not required to pass precisely 
through the data values. For the former process errors 

in t(p) are of the same order as the errors in the T(X) 


values. Kennett (1976) found that values of t(p) obtained 
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Figure 3.3 


Travel-time curve of first branch from the 


Yukon line of Project Early Rise. 
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Figure 3.4 


The function 1T(X,p) for the first travel-time 
branch from the Yukon line. of Project Early Rise.) The 


p values are such that Psa ase 
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from the latter process were much more widely scattered 
than estimates using the graphical construction. Spline 
smoothing tends to overfit the data to a degree which is 
dependent upon the nature of the spline routine and I 
recommend an alternate procedure in which error estimates 
result directly from the process. 

For real seismic data, travel-time branches are 


Suitably approximated by a polynomial of the form: 


Ts gq +a bDXK+ oe 2 (3.24) 


For a branch consisting of N points (T,,X5, i=1,N) the 
quality of the parameterization (3.24) is given by the 


root mean square error, 


N 1. 
2 2 g 
E(a,b,c) = ie (adpby, tock. = Ty Vs (35258 
1i=1 
Prom (37)).,. (3.24), and the condition that ay = 0 the 


expression for the branch in the t(p) plane is 


2 2 
he Ae ag: 1 ee 
Ap Pee wee Oat ei 
b + 2xC, EA Ope es b + 2cX, ;: Gyx 2 (3..26)) 
Dig ACK sep apt Jenne c > 0 


where Xy and xy are the smallest and largest X coordinates 


of the branch respectively. The best fit solution for 
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a= aor b= Dor CG for which E(a,b,c) is a minimum. 

‘ Gt canals). pola, , 
Using aa? ecm 0 and (3.26) the best solution is 
defined by 

ay Ay A. A3 By 
De = |A, A, Ay Bo 
igs gee epoa waleagelsy rushes 
2 
N i.e Xs, pnd hae 
i i i 
2 3 
=a ie nD a DG Be, ae ioe (3.27) 
i i ods 
cle ee Be Ex?T, 
x a i ree 
with solutions taken over the points i = 1,N on the branch. 
Now, letting Eo = E(ajrborc,) and Ey = SEor consi- 


der the set of real values (a,b,c) such that E(a,b,c) s Eye 
Clearly when €&<1, the set is empty, and when € = 1, the 
set consists of one point in (a,b,c) space, namely a = aoe 
b = Doe aig oie For the case of € > 1, the set consists 
of more than one value (a,b,c) and it is possible to show 
that the set is of full measure. This multi-solution 
nature of the parameterization of T(X) data has in past 
caused ambiguity in the interpretation of seismic obser- 
vations. Results from many long range refraction experi- 
ments have been interpreted in terms of layers of constant 


velocity; travel-time branches are parameterized by the 


simpler form T = a + bX, and the layered model is inferred 
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from the values of "a" and "b" from the various branches. 
Clearly, small errors in the values of "a" and "b" 
change the deuce layered model making comparison of 
velocity models from different studies difficult. 

There is no need for confusion however, since it 
1S possible to find all values (a,b,c) which fit the T(X) 
data within any root mean square error E, > E_. The 


e Oo 


solution space corresponding to E, > Ey is composed of 


. 
all points within the ellipsoid defined by, 


2 2 
A,a 4 A,b “f 2A,ab + 2a (cA, - B,) + 2b (cA, - B.) 
Pea. Senta Abe aCe ein = 20 (3.28) 
5 3 6 O 1 : 
where Ag 2 BT sh The ellipsoid is bounded by the planes 
= git and c = Coase 
Bey gs ot HO endo) 21/120; ) 
Cmin ~ ~ \22 29 2423 2] 
(3.29) 
Soaks) (0. Mote 1/126.) 
Cmax ~ 2 2 123 21 
where 
Geena A ge oan Bl ARN Ry ok oR nh at Re 
he eilinee” 4 = 3°°4 2°°5 °3 
« 2 2 
0, = 2(8,(,-454,) +8, (A, A, > ALDA,) + Bo (AS = ALAS) 


e 2 2 2 | 
QO, = (Ag- (E,) A,) (A,A,- A) +B, (2A,B,- A,A))- BZA, . 


THYS BrOVLded 6... «kc UX. 1c , the solution space in the 
min max 


(a,b) plane are points contained within the ellipse given 
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ats 223 ),. 

In general for a fixed value of c, the ellipse in 
the (a,b) plane is of the form shown in figure 3.5. For 
typical experimental data the major axis is tilted clock- 
wise from the "a" axis by a few degrees and the major 
axis is much greater than the minor axis (the ratio of 
the major to the minor axis is much greater than that 
shown by figure 3.5). The reason for this ‘tilt’ can 
be understood Be considering the fact that experimental 
T(X) values for a given branch fall almost on a straight 
line with intercept ay and slope bo: If an attempt is 
made to fit the data to another straight line for which 
the time intercept a 7 ae is fixed, the resulting value 
of the slope, b, will be such that b < bos Similarly if 
saan then b > bo. Also, the ellipsoid defined by (3.28) 
is such that when c increases the allowable values of b 
decrease. These properties of the solution space are not 
accidental; they are related to the stability of the 
function t(p). For example, consider figure 3.5 and 
assume; for simplicity, that c= 0. For-any (a,b) within 
the ellipse we have from the mean square error in (3.25) 
that t(p) = a where p = b. Because of the tilt of the 
major axis of the ellipse, points with higher b values 
have lower a values; this situation is in accordance with 
the monotonic decreasing behaviour of tT(p). Also, since 


the ellipsoid encompasses smaller b values for larger c 
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Typical error ellipse in the (a,b) plane for 
T = a+ bx+ cx? type parameterization of travel-time 
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values, the tendency is toward preservation of the range 
of p defined by (3.25) when c changes. 

In order to find all possible expressions of the 
T(X) branch in the t(p) plane for a given error BE or we 
need only calculate t(p) from (3.25) for points (a,b,c) 
on the surface of the ellipsoid defined by (3.28). The 
image curve t(p) of any interior point (a,b,c) will be 
"contained! by the image curves of the surface points. 
This fact can be seen by aera acting the dashed line in 
figure 3.5; from (3.25) the images t(p) of points on the 
- dashed tee have the same p range and if T, (p) and To (Pp) 
result from the two endpoints then T, (p) <T34 (P) <To (Pp) 
where Tse (P) is the image of any point between the end- 
points. The procedure then is to calculate t(p) from 
(3.24) for discrete points on the boundary of each 
ellipse corresponding to discrete values of c between 
Caan and Canes With the solution space defined by (3.28) 
we can calculate all possible t(p) expressions of every 
T(X) branch for any root mean square error. An envelope 
in the t(p) plane can be determined by examining all 
possible solutions. Finally this envelope may be mapped 


into an envelope in the V(Z) plane by using the technique 


described by Bessonova et al (1974). 
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Results from Project Early Rise 


The travel-time branches recorded in Western 
Canada from shots in Lake Superior along the 'Yukon' 
line during Project Early Rise together with the result- 
ing t(p) envelope are shown in figure 3.6. The p axis 
of the t(p) plane has been normalized to a surface velocity 
of 4 km/sec. The letters indicate the location of branches 
in the t(p) plane. Each travel-time branch was fitted to 
a best fit curve of the form given by (3.24); the uncer- 
tainty ellipsoid (equation 3.28) in (a,b,c) space for 
each branch was then determined for a value of r= ee 
The upper and lower t(p) bounds for each branch were 
determined by sampling the surface of the error ellipsoid; 
the p range was restricted to the region defined by the 
best fit curve (3.24). With the t(p) bounds determined 
Eom each (branch an this manner it was found that the 1t(p) 
regions for various pairs of travel-time branches blend 
into one another; thus the t(p) images of branches C and 
D have a small common region along the p axis, as do the 
images from the pair E,F and the pair G,H. The pieces 
of t(p) envelope with common regions mentioned blend into 
one another quite smoothly and therefore it was not 
difficult to form the 'composite' t(p) upper and lower 
bounds in the intersecting p regions. This phenomenon 
of blending in the t(p) plane is not surprising in view 


of the fact that the travel-time branch pairs involved 
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Figure SO 


The insert on the upper right-hand corner shows 
the travel-time branches for the Yukon line recorded 
during Project Early Rise. The main diagram shows the 
Tau envelope for this data. The letters indicate loca- 


tion of branches in the tT(p) plane. 
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do in fact converge in the T-X plane. Now there is also 
a region in the t(p) plane for which the p regions of 

two branches overlap but for which the values of 1t for 
the overlapping p interval are offset; this jump in 1t(p) 
is indicative of a low velocity zone. The t(p) envelopes 
resulting from branches D and E have a common p region 
defined by 0.465 < p < 0.478. The maximum jump in 1T(p) 
within this interval is defined by the maximum of the 
difference between the upper 1t(p) bound for branch E and 
the lower t(p) bound for branch D; the value is appro- 
ximately 1.9 seconds. Similarly the minimum jump in 1T(p) 
in this p interval (the minimum of the difference between 
the lower t(p) bound for branch E and the upper t(p) bound 
for branch D) iS approximately 1.1 seconds. Thus the 
parameters which define the low velocity zone kinemati- 
cally are q, = 0.465, qa = 0.478, og, = 1.1 sec, and 

oy = 1.9 sec. Now the t(p) envelope for p regions 
between the images of the individual branches has been 
constructed using the fact that t(p) is a monotonic 
decreasing function. Notice that the t(p) envelope in 
these gap regions is responsible for the major uncertainty 
(that is the t(p) envelope is widest for p intervals 
prasedit branches A and B, B and C, and F and G). Note 
also that, of the rays observed, those having the least 


depth of penetration constitute branch A. Thus T(p) 


data 'commences' at ap value of about 0.61 and there 
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are no ‘ovservations Lor the,range 0.64. <p <) 1.0, The 
extent of this interval may seem overwhelming but it 

is rather artificial in view of the assignment of the 
generous value (4.0 km/sec) for the surface velocity. 
The t(p) envelope in this region has been constructed 

by linear interpolation between the existing point 

Ty (0.61) and 1(1.0) = 0 for the upper bound and linear 
extrapolation between the existing point Tt, (0.61) and 
the point 1(0.66) = 0 (note that the value of p = 0.66 
corresponds to a surface velocity of 6 km/sec) for the 
lower bound. It has been mentioned that the p ee of 
each branch in the t(p) plane has been dictated by the 
values of b= bo and c = om (that is the coefficients 
describing the best fit curve of the form (3.24)). It 
is worth noting that if any other values of "b* and ‘c’ 
within the error ellipsoids are taken then the resulting 
image t(p) curves essentially fall within the T(p) 
envelope shown in figure 3.6; any points which fall 
outside the t(p) bounds shown are so close to the boundary 
that hey are. 1nslgnificant. 

The t(p) envelope in figure 3.6 has been inverted 
using the Bessonova et al (1974) method described in the 
preceding section. A low velocity zone with parameters 
dye dye d,, and os specified above has been assumed. The 
resulting velocity-depth envelope is shown in figures 


3.7(a) and 3.7(b). The darkened region in figure 3.7 (a) 
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Figure 3.7 


Galooity bounds for the Yukon line t(p) bounds. 
In part (a) the darkened region is a barrier such that 
no velocity-depth curve may pass entirely through it 
without first exhibiting a low velocity zone. The 
dashed lines in part (b) outline possible low velocity 
zone configurations when the top of the low velocity 


zone is assumed to be at a depth of 60 km. 


VELOCITY (KM/SEC) 


VELGCITY CKM/SEC) 


8 


6 


200 400 
DEPTH (KM) 
200 400 


DEPTH (KM) 


600 


600 


gs S 


ee sa hhh, q 
7 ore i ne e ai 
a uD : 
aiienies 7 A ay u 


a ) von ‘ 
7 
. ye ae 
1, CA ie at Oe 
ne i 7 aM ¢ 
7 oe, att 
‘ie a vy 
: . \ 


wittamot first vt sxhinceing! a tow 4 
dsulie inefdin part tb) outline & 8 


yee atta eetarhaun when the top 


i] 


ia" 


\ / ‘, 


Cibo ara Naan 


(8) 


156 


is a 'barrier' associated with the low velocity zone. 

The darkened region has the following properties: 

(a4) No velocity-depth curve may pass entirely through 
it without first describing a low velocity zone. 

(ii) The low velocity zone must commence at some 
‘valid’ point which is contained within the darkened 
region. The low velocity zone has a magnitude defined 

by o, where Cake: oy. Now in the inversion it was 
implicitly assumed that the maximum (flat earth) slowness 
attained in the low velocity zone was p = 1.0. This 
slowness corresponds to a spherical earth velocity of 
approximately 4.0 km/sec for the depth shown. In view 
of this minimum bound (4.0 km/sec) for the velocity 
within the waveguide it is clear that for any allowable 
value of o there will be points within the darkened 
region which are so deep that it is impossible to  cons- 
truct the corresponding waveguide commencing at these 
points. Any point within the darkened region from which 
it is possible to construct the low velocity zone for 

BM Jy Seago Lite Odin + 

Figure 3.7(a). shows that the top of the. low velocity zone 
may be as shallow as 60 km, and also that the waveguide 
may be as thick as about 150 km (in this case the 
difference between the velocity within the waveguide 

and the velocity at the top of the waveguide is infinite- 


simal). Figure 3.7(b) illustrates the low velocity zone 
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configuration for o = o, when the top of the waveguide 


1 
is at a depth of 60 km and the velocity at the top of 

the waveguide is about 8.3 km/sec. The differences in 
depth between the upper and lower dashed lines is appro- 
ximately the thickness of the low velocity zone as a 
function of the velocity within the zone. This rela- 
tiensiip 1S Crue. if Tt is assumed that. the “ftlati earth ' 
waveguide is rectangular in which case the corresponding 
spherical earth waveguide (equation 3.2) is approximately 
rectangular. The low velocity zone configuration for 

= oF Ts very Close, to thae for vo, = O45 and thus dashed 
curves similar to those in figure 3.7(b) for the case 

Oo = 9, are not shown. If we consider a velocity model 
which lies midway between the upper and lower bounds 

Shown in figure 3.7 then the 'form' of the bounds suggests 
a crustal thickness of about 30 km, a rapid increase in 
velocity near about 170 km depth (as the velocity curve 
exits the low velocity zone), and another rapid increase 
in velocity near: about 420 km. The presence of the upper 
mantle structures is not a new phenomena. Indeed continen- 
tal P wave studies which reveal the asthenospheric low 
velocity layer, a high velocity gradient immediately 
beneath the low velocity zone and a high velocity gradient 
near about 420 km depth include those of Lukk and Nersesov 


(1965), Niazi and Anderson (1965), Johnson (1967), Julian 


and Anderson (1968), Green and Hales (1968), Archambeau 
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et al (1969), Helmberger and Wiggins (1971), Hales (1972), 
McMechan and Wiggins (1972), Wiggins and Helmberger 
(1973), Massé (1973), and Massé and Alexander (1974). 

In addition the P wave travel time data of Dowling and 
Nuttli (1964) and Hill (1972) indicated an asthenospheric 
low velocity zone; high velocity gradients near 420 km 
were reported by Golenetski and Medvedeva (1965), Barr 
(1967), Lewis and Meyer (1968) and Mereu et al (1974) on 
the basis of continental P wave travel times, If is of 
particular importance that this study has revealed a well 
developed low velocity zone associated with a stable 
platform region. The low velocity zones detected in the 
majority of the studies cited above are associated with 
tectonically active regions. Some studies, however, 
suggest that associated with the asthenosphere underlying 
the stable central and eastern United States there is a P 
wave low velocity zone which is albeit not as pronounced 
as the one underlying the tectonically active western 
United States. Thus Hales (1972) cites the presence of 

a low velocity zone of large magnitude in the western 
United States and one of lesser magnitude in the central 
and eastern United States. Also Dowling and Nuttli (1964) 
find differences between the low velocity zone structures 
in the western and eastern United States; they state that 
the top of the low velocity zone is at a depth of about 


70 km in the west and about 125 km in the east. Massé 
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(1973) studied travel time data from various lines of 
Progect Han lyeriserend Tounda low velocity. 2Zonesin) tne 
depth range 94 to 107 km with an average velocity of 
8 km/sec (which is about 0.4 km/sec smaller than the 
velocity at the top of the low velocity channel) in the 
central and eastern United States. The parameters (depth 
and magnitude) of the low velocity zone described by 
Massé (1973) agree very favourably with results from 
this study (figure 3.7). Furthermore Massé and Alexander 
(1974) concluded that the upper mantle P velocity beneath 
Scandinavia and Western Russia closely resembles the model 
derived by Massé (1973) for the central and eastern United 
States. Thus this study and other studies have revealed 
the presence of substantial low velocity zones associated 
with the asthenosphere underlying stable platform type 
regions. This fact lends credence to the prediction and 
tenets of Hales (1972) who states "... I hazard the guess 
that further study will show that the low-velocity eone 
is a general and required feature of the upper mantle... 
ifs) quite, clear that tthe: relative motions of the: outer 
part of the earth, which may now be regarded as established 
beyond reasonable doubt, do not take place only in tec- 
tonicallys active areas," 

Now it is very interesting to compare the results 
of this experiment with the results of various other re- 


fraction surveys conducted in North America. A comparison 
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of results in the velocity-depth plane gould be mislead- 
‘ing since in past various authors have used a variety of 
inversion techniques. It is felt that the tT(p) plane 
offers the most unbiased te for comparison or the 
kinematic properties of crustal and upper mantle P waves 
which eraverde various regions of the North American 
continent. Iyer et al (1969) list the travel times and 
distances for travel-time branches from Project Early Rise 
refraction lines. For the purpose of enegesrm with this 
experiment, these travel-time branches ee fitted to best 
fit curves of the form (3.24); the corresponding best fit 
t(p) Curves were then determined using equation (3.26). 
The lines used and fhe approximate distance ranges of the 
travel-time ee eee considered are as follows: Chapleau, 
Ontario to Chibougamau, Quebec (510 to 1080 km), 
Chibougamau, Quebec to Glace Bay, Nova Scotia (1100 to 
2040 km), Saint Cloud, Minnesota to Holyoke, Colorado 

(440 to 1330 km), Bemidji, Minnesota to Miles City, Montana 
(370 to 1280 km), Dyment, Ontario to Sheho, LEV eC ESE 
(620. tO BLIEO km), Deduc, Alberta to Post. 140, ee alert 
Columbia (2150 to 2490 km), Fort Nelson to.Post 750, 

By Leash Columbia (2540 to 2680 km), Flin Flon, Manitoba 

to Yellowknife, Northwest Territories (1380 to 2310 km), 
Bear Lake, Michigan to Potters Hill, North Carolina (40 to 
1680 km), Lake Superior, Michigan to Eldorado, peeisae 


(440 to 1430 km), Lake Superior, Michigan to Cornstock, 
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texas: )(l/0) tor2250 km), Colorado Springs, Colorado co 
Charlston, West Virginia (1110 to 1620 km), Chapleau, 
Ontario to Chibougamau, Quebec to Schefferville, Quebec 
(460 “to 1700 km), South Eastern Ontario (300° to*7 00. km)), 
Lake Nipigon to Smooth Rock Falls, Ontario (240 to 610 km), 
Port Arthur, Ontario to Baralzon Lake, Manitoba (300 to 
1490 km), Bemidji, Minnesota to Vancouver, British 
Columbia (470 to 2600 km), and from Lake Superior to the 
LRSM stations (500 to 3480 km). Thus the data set sam- 
ples central and eastern North Américan regionsgand the 
distance ranges include those of this experiment. Points 
on the best fit t(p) curves for the travel-time data 
mentioned above are shown in figure 3.8 together with the 
T(p) envelope derived from the 'Yukon' line of this experi- 
ment. The agreement between the t(p) points from the 
various lines and the t(p) envelope shown is quite remark- 
able; with the slowness scale of figure 3.8 adjusted with 
p = 1 corresponding to an apparent velocity of 4.0 km/sec 
it can be seen that the maximum departures in apparent 
velocity of points from the envelope boundaries are about 
0.4 km/sec near 8 km/sec, 0.2 km/sec near 8.4 bn/aee and 
about 0.2 km/sec near 9.6 km/sec. 

Now let us consider North American seismic refrac- 
tion data which have been interpreted by various authors; 
in, particular we Shall again examine results in the t(p) 


plane. Not all published reports list data from which 
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Figure 3.8 


The Yukon line t(p) envelope with t(p) points 
based upon the Project Early Rise travel times listed 


‘by Iyer et al (1969). 
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t(p) information can be conveniently extracted but it is 
felt that the following studies taken collectively pro- 
vide a reliable overview of interpreted North American 
seismic refraction surveys; the authors and the geo- 
graphical areas studied are Dowling and Nuttli (1964)- 
eastern and western United States, Roller and Jackson 
(1966) - central United States, Hobson et al (1967) - 
Hudson Bay, Barr (1967) - Canadian Shield, Ruffman and 
Keen (1967) - Hudson Bay, Green and Hales (1968) =- central 
United States, Lewis and Meyer (1968) - central United 
States, Gurbuz (1970) - central Canada, Mereu and Jobidon - 
eastern Canada, Hill (1972) - western United States, 
McMechan and Wiggins (1972) - western United States, 

Berry and Fuchs (1973) - north eastern canadian Shield, 
Massé (1973) - central and eastern North America, Wiggins 
and Helmberger (1973) -~ western United States, and Massé 
(1974) - central and eastern North America. Points in the 
t(p) plane from these studies are shown collectively in 
figure 3.9 together with the t(p) envelope from this 
experiment. Unlike the t(p) points due to the travel- 
time listings of Iyer et al (1969) (figure 3.8), there 

are several points which diverge considerably from the - 
t(p) envelope in figure 3.9. There are two major reasons 
for this disparity. Firstly the physically ee ehes iat 
negative values of t(p) in the p range of about 0.62 to 


0.67 (apparent velocities from about 6.0 to 6.5 km/sec) 
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Figure 3.9 


The Yukon line tT(p) bounds and t(p) points 
based upon the results of many seismic surveys 


conducted in North America reported by various authors. 
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are a result of the Hudson Bay observations of Hobson 

et al (1967) and Ruffman and Keen (1967). The combined 
effects of a lateral inhomogeneous Hudson Bay crustal 
structure and inconsistent navigation can cause this dis- 
crepancy (indeed Ruffman and Keen (1967) examined the 
consistency of three possible navigation series). The 
second and more interesting reason for the disparity is 
that the data set shown in figure 3.9 is composite, being 
due not only to surveys conducted in central and eastern 
North America but also to studies applicable to the 
tectonically active western United States regions. 
Anomalously large t(p) values ("anomalous' if compared 

to the t(p) envelope) near p = 0.61 (apparent velocity = 
6.6 km/sec) and in the p range of about 0.46 to 0.51 
(apparent velocities of about 7.8 to 8.7 km/sec) are due 
to the western United States investigations of McMechan 
and Wiggins (1972), Wiggins and Helmberger (1973), Hill 
(1972), and data from the western United States line 
studied by Dowling and Nuttli (1964) (Dowling and Nuttli 
(1964) also discussed results from the central, eastern, 
and southern United States and the t(p) data from these 
areas fall within the t(p) envelope or close to the 
boundary). The difference in kinematic properties between 
rays which traverse the tectonically active western United 
States and those which traverse the most stable regions 


of central and eastern North American regions is not 
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Surprising. Although a determination of an exact velocity- 
depth function requires complete specification of T(p) 
it can be stated that the high t(p) values for the 
western United States in the p range of about 0.46 to. 
0.51 (figure 3.9) are consistent with 'overall' lower 
velocities in that tectonically active region as com- 
pared with the upper mantle velocity of the central and 
eastern portions of North America; thus rays which pro- 
pagate through a higher magnitude low velocity zone and 
bottom just below it will suffer larger delay times and 
thus larger t(p) will result. Other points which lie 
far from the t(p) envelope in figure 3.9 are located near 
p = 0.47 and t = 4 sec; here the low p values as com- 
pared with the t(p) envelope are due to the study of 
Berry and Fuchs (1973) (some of their points also fall 
within the t(p) bounds). The spread in p values of the 
T(p). duestoS%Berry and“ FuchsyY (1973)“as not surprising 
Since the survey area included various lines crosSing a 
large gravity anomaly associated with the Superior- 
Churchill province boundary in the northeastern Canadian 
Shield. Finally it is emphasized that the t(p) points 
in figures 3.8 and 3.9 which are applicable to central 
and eastern North America fall within the t(p) envelope 
or very close to the boundaries. Also deviations of the 
t(p) points from the envelope and the spread in p values 


of “the t(p) points in’ figures 3.8 and. 3.9 become less 
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pronounced as p decreases indicating that the upper 
mantle structures of the various regions sampled tend 
to coalesce with aes Figures 3.8 and 3.9 also show 
t(p) points for rays which have bottomed at greater 
depths than those of this study. It is not intended 
that the comparisons discussed here constitute an 
exhaustive analysis of all refraction surveys ever 
conducted in North America but rather it is hoped that 
the discussion illustrates the importance of the tT (p) 
plane. 

A final point regarding the width of the velocity- 
-depth envelope in figure 3.7 is worth mentioning. It 
can be seen from figure 3.6 that the slowness extent of 
the t(p) envelope for t(p) values from about 12 to 19 
seconds is about Ap = 0.13 which corresponds to an uncer- 
tainty in apparent velocity of about 0.4 km/sec. On 
the other hand the velocity envelope in figure 3.7 
defines a velocity uncertainty of about 0.6 km/sec in 
the corresponding depth range of about 200 to 360 km. 
The velocity-depth boundaries shown in figure 3.7 define 
a rather large velocity uncertainty partly because the 
corresponding t(p) envelope bounds a finite area (we do 
not know t(p) exactly) and partly because the inversion 
scheme described in preceding sections does not yield 


greatest upper and least lower velocity depth bounds 


corresponding to given t(p) bounds. In view of the nature 
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of travel-time observations (observations are incomplete 
and in particular it is difficult to extract seismic 
signals which correspond to rays bottoming in regions 

of low velocity eae from the seismogram) the inver- 
sion technique described cannot yield detailed velocity 
depth information. It is hoped that in future the exten- 
sion of the Tau method to a Tau inversion method which 
incorporates partially reflected and converted phase will 
tend to alleviate this problem. The consideration of 
amplitude characteristics should also result in narrower 
velocity depth bounds; for example,McMechan and Wiggins 
(1973) used the Cagniard-de Hoop technique in conjunction 
witha Hexglotz-Wiechert-type.extremail , inversionascheme 
to construct models of the P wave structure for the 


upper mantle. Also Bessonova et al (1976) have described 


a Tau inversion method which yields more accurate velocity 


depth bounds than the one reported by Bessonova et al 
(1974). There is a further complication however ~- and 


that is the effect of lateral inhomogeneities on travel 


times, amplitudes, and the vector ray parameter in studies 


of the upper and lower mantle and the core of the earth. 
Possibly future two dimension array studies in conjunc- 
tion with the above mentioned techniques will increase 


our knowledge of the detailed structure of the earth. 
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Conclusions 


Not only is the Tau method an objective approach 
to the SPonten of inversion of real seismic T(X) data 
sar cape we it may be used to examine the variation of 
Groertalnty of the velocity depth function due to changes 
in seismic array geometries. In the case of real seismic 
data, estimation of the function t(p) from T(X) can be 
performed satisfactorily using the 'Clairaut' approach 
as long as errors in T(X) are inconsequential. With 
realistic scatter, due to inhomogeneities or errors in 
measurement as is the case for long range refraction work, 
it becomes necessary to smooth the T(X) data in order to 
derive a meaningful velocity structure. An unbiased 
method of smoothing involves a determination of all 
possible curves, T = a + bX + cx’, that fit teach: T(x) 
branch within any prescribed root mean square error. The 
solutions (a,b,c) are found to be contained within an 
ellipsoid in (a,b,c) space. ‘The properties of ‘the ellip- 
soid, for real data, reflect the inherent stability of 
the function t(p). 

Application of the Tau method of inversion to 
travel-time data applicable to western Canada indicates 
a crustal and upper mantle P wave velocity structure 
which consists of a crust approximately 30 km thick, a 
low velocity zone between depths of about 60 and 210 km, 


and rapid increases in velocity near 170 and 420 km. 
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A comparison of the t(p) bounds derived for western Canada 
with t(p) points associated with other regions of North 
America reflect crustal and upper mantle lateral inhomo- 


geneities. 
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Appendix 1 


A vector norm is a meaSure of the size, in some 
Sense, (Ob a Vector. Live is a vector then the most common 
vector norm is the Euclidean norm, giving the length of 


the vector, 
Nes exe: 


Any vector norm has the following properties: 


(i) pics| T= Meter xa Mtoreald urea lc. 

(iia) Dichie Of ste xt is! norammiLlavector: 

Serr ehy ies xl it live cone al lavectons a.and yA), 
The induced matrix norm, ||A|| of ann x on matrix, 


A, is given by 
[Jal] = max] [ax] | where Pear resales 


In this equation x has been normalized and for the 
Euclidean vector norm, ||A|| is called the ‘spectral’ 
matrix norm and is equal to the length of the longest 
vector in the image set {Ax} under the transformation 
Koy = Ase. Matrix norms satisfy the properties (i) to 
(iii) above. It is possible to show that the spectral 
norm |jal| is equal to the maximum singular value of A 
(equivalently Reels is equal to the spectral radius, 
DRE G Of the matrix ata which is the non-negative 


Square root of the necessarily non-negative largest 
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eigenvalue of the matrix AC nye Another useful vector 


norm is the 'norm-infinity' vector norm defined by 
eles) | eee eax! a fc 


where 


e T 
ES NES geen aise 1S 


The induced norm-infinity matrix norm is 


N 
[IA] |, = max eee | 
1<i<N J=1 <2 
where 
A= eae : bigs ) Sesh ye hk, ik Ome 
The ie daa on number', K(A), of square non-Singular 
matrix A is defined by K(A) = }jal].|Jar?] PG leae hy i.e 


condition number depends upon the matrix norm which is 


chosen. Regardless of the matrix norm, however, it can 


Be@shown that. [|b ("Se (a) P< 2Sheror the™@spectral norm’ « (A) 


Up Ue where Uy and UW, are respectively the largest and 


smallest, necessarily positive, eigenvalues of the matrix 


lb 


AN 2° "in ¥géneral*k(A)-*is avmeasure "of the "stability ‘of the 
inverse problem given by equation (1.2). Unstable systems 


are associated with large condition numbers and are termed 


'i311-conditioned'. 
A final point worth mentioning is that if y = Ax 


then for any vector norm and associated induced matrix 
norm, 
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Appendix 2 


The analysis of diagonal matrix form under ortho- 
gonal equivalence is based upon the following theorem and 
GUscussion..« If ‘Avis’ a real -nixn matrix then there exist 
two real@nex<n forthogonal-matrices, U,V so that utav isa 
diagonal matrix D; furthermore we may choose U and V so 


that the diagonal elements of D are 
HW, 2 U52 bd Sl Ta 2 eed eee = u = 0 , 


where r is the rank of A. The numbers Uzpreeer, are the 
Singular values of A. Now consider A as representing a 
linear transformation of one n-dimensional space G into 

a second such space H. Thus h = Ag is in H for any g in 
G. Now effect an orthogonal change of coordinates in 
Space G and a different orthogonal change in space H such 
that we have the new representations ts Vg' and h=Uh) 3 
As a result of the changes of bases in G and H the trans- 
formation A obtains the simple representation of the 
matrix D. Thus we have h! = Dg'. The transformation 

now Simply maps the first coordinate axis of G onto the 
first, coordinate axis’ of H with a pear ton ime lenetay 

Uy > 0. It does the same for the ee SES, att ee coor- 
dinate axes with magnification factors Horeeer HL. The 
(r+1)-th,...,n-th coordinate axes of G are mapped onto 


the zero vector of H. The most important result for the 


purpose of applications to seismic arrays is that D maps 
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the unit hypersphere Ss = {g':||g'||=1} into an r- 
dimensional hyperellipsoid E = DS of vectors h' such 


that 
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Appendix 3 


Tau Inversion Program 


The following computer program calculates velocity- 
depth limits for given t(p) limits using the algorithm 
given by Bessonova et al (1974). It is assumed that 
the t(p) limits are calculated from the normalized travel 


time, I, and distance \observatnons, x", “that. is 
BL, =A / 2 


et 


X/2 VSUR 


where T(sec) and X(km) are the actual travel times and 
distances and VSUR(km/sec) is the velocity at the surface 


of the earth. Thus 


B{D)asjT> oepx' 
and 


Oe Fad Jax 


The upper and lower bounds Ty (Pp) and Ty (p) fox the func- 
tion t(p) are given in terms of piecewise second order 


polynomials in p and 


MU 


II 


number of polynomials defining Ty (P) 
ML = number of polynomials defining Ty (p) : 


For the function Ty (P), the p axis is segmented by the 


values 
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NU =4 MU 41, PUG) a= « PO 


and 


PU(I+1) < PU(I) LOLA Ee 


The function Ty (P) is then defined by, 


Ty(p) = AU(I) + BU(I)p + CU(I)p* for PU(I+1) < p < PU(I). 


Similarly for the function tT, (p) the p axis is segmented 


by the values, 


POLIS) £ St1s NE 
where 

NL = ML +1, PEL) = 0 
and 


PRG) vei) stor -abio a3 


The function tT, (p) is then defined by 
1, (p) = AL(I) + BL(I)p + CL (I) p* for PICT <a = BL or 


It is assumed that PU(NU) = PL(NL). Now NLVZS is equal 
to the number of low velocity zones under consideration, 
The quantities QU(J) and QL(J), J = 1,...,NLVZS define 
the region of uncertainty along the p axis for which T(p) 
may exhibit a 'jump' corresponding to the Jth low 
velocity zone. Thus for the Jth low velocity zone the 


jump may take place at a slowness value q for which 
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QL(J) < q < QU(d) 


It is assumed that 


OO 40m) ) 
and 


OUT Jee or 0 : 


The quantities SIGU(J) and SIGL(J) define the uncertainty 
of the jump in t(p) associated with the Jth low velocity 
zone. Thus if At is the jump in t(p) associated with 


the Jth low velocity zone then, 


SRG yes it So Laue Ss 


The quantity UMAX(J) is the maximum slowness which can 
be attained within the Jth low velocity zone. 

The values SAM and NDIV define coarse and fine sampling 
parameters respectively used in the search for upper 

and lower velocity bounds. SAM is a slowness sampling 
used to determine a rough estimate of the location along 
the p axis for which upper and lower bounds are attained; 
a more precise estimate is then attained with a search 
conducted in the region of the rough estimate with a 


slowness sampling interval of SAM/NDIV. 


Input 


The input must be in the following order and 


Format. 
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NU, NL FORMAT (215) 
PUACTOe iL lee OND FORMAT (2F15.3) 


BiG 7 ok 


i 
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FORMAT (2F15.3) 


(ALL(L) BUG) pCUKT),) ¢ I= FORMAT (3F20.0) 


| 
ke 
= 
a 


CART), BE (HE) 7eCR (1) ) Yur FORMAT (3F20.0) 


| 
ke 
= 
BR 


Pi) Sd MU FORMAT (1015) 
If CU(I) = 0 then ILU(I) = 1, and otherwise 


LLU GT ips “Os 


ie (EPs ar = Te FORMAT (1015) 
If CL(I) = 0 then ILL(I) = 1, otherwise 


PL, Ce he 0 - 
SAM FORMAT (F15.3) 
NDIV FORMAT (15) 


VSUR FORMAT (F15.3) 
VSUR is the actual velocity at the surface of 


the earth in km/sec. 


NREGS FORMAT .(1i5) 
NREGS is related to the number of low velocity 
zones; NLVZS . If QL(NLVZS) > PU(NU) then NREGS = 
2 NLVZS + 1. If QL(NLVZS) = PU(NU) then NREGS = 


2 NLVZS,. 


NVELS (K), K = 1, NREGS FORMAT (1015) 
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This quantity gives the number of slowness values 
for which upper and lower depth bounds are acquired 
within certain regions of the slowness axis. For 
this purpose the lst region of the slowness axis is 
mniveny by "OU (1) “<°p-s 1. “The Kth region for K > 1 and 
K even is given by QL(K/2) < p < QU(K/2). The Kth 


region for K > 1 and K odd is given ‘by 


DUAKED) / 2). 08 ys OL C(Kel )/ Zac al as 


If there are no low velocity zones then there is only 
one region, namely PU(NU) < p < 1. Upper and lower 


depth bounds are calculated for equi-spaced slowness 


values within the regions. 
13. NLVZS FORMAT (15) 


14 SsuiGHCl (de, SEGH W) ,0U (5) AOL (a). DMA) MUSTAGH eT = pNLVZS . 


FORMAT (5F15.5,15) 


The previously undefined quantity MUST(J) refers to 


the Jth low velocity zone. If MUST(J) = 1 then a low 
velocity zone must occur and if MUST(J) = 0 then a 
low velocity zone may occur. Notice that MUST(J) = 0 
Th ait ORLY “Lo orci) 80 z 

Output 


For each low velocity zone the following quanti- 


ties are given, 
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(a) the maximum slowness of the zone assuming the 


maximum thickness and t(p) jumps of SIGL(J) and SIGU(J). 


(a>) the maximum thickness of the zone assuming the 
maximum possible slowness and t(p) jumps of SIGL(J) 


and SIGU(J). 


The above results are given in terms of the 'flat earth" 
velocity depth bounds associated with the normalized 
tT(p) envelope. 

Upper and lower depth bounds for velocity are 
given in terms of an actual spherical earth; depths are 


given in km and velocities in km/sec. 


Operational Procedure 


The calculation of the upper depth limit is done 


in SUBROUTINE UPPLIM which in turn calls SUBROUTINE ULIM. 


Similarly lower limit calculations are performed by 
SUBROUTINE LOWLIM which calls SUBROUTINE LLIM. Low 
velocity zone calculations are done within the main 
program. At present the program is set up such that it 
can handle up 5 low velocity zones, 209 piecewise poly- 
nomial descriptions of the upper and lower t(p) bounds, 
and give 500 estimates of upper and lower depth bounds. 
Tf more of the above parameters are required then a 
simple change of the dimension statements would be 


necessary. 
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31 


38 


39 


41 
42 


43 
100 
101 
102 
2222 
2223 
4343 


4345 
4346 
7000 


2228 
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DOUBLE PRECISION SIGU(5),SIGL(5) ,QU (5) ,QL (5) 
DOUBLE PRECISION HMAX (5) ,HMIN(5) , UMAX (5) ,UMIN(5) 
DOUBLE PRECISION A(500),B(500) , ABEST (500) , BBEST (500) 
DOUBLE PRECISION VA (500) ,VB(500),SAVE (5,50) 
DIMENSION MUST(5) ,NVELS (11) 

DOUBLE PRECISION SKE, TH, THM,SLOW,PD, ARG1,ARG2 

DOUBLE PRECISION PU (210),PL(210) 

DOUBLE PRECISION AU (210), BU (210) ,CU (210) 

DOUBLE PRECISION AL (210) , BL (210) ,CL (210) 
DIMENSION ILU (210) ,ILL(210) 

DOUBLE PRECISION DIV, SAM, VEL, DEP, R, VSUR, FACT 

COMMON PL, PU 

COMMON AU,BU,CU, AL, BL,CL 
COMMON SAM 
COMMON NDIV,NL,NU 
COMMON ILL, ILU 
FORMAT (1H ,*THE MINIMUM POSSIBLE SLOWNESS IS? 
1, F15.10) 

FORMAT (1H , 

1" MAX SLOWNESS FOR MAX THICKNESS AND MIN SIGL=' 
1, F15. 10) 

FORMAT (1H , 
1*MAX THICKNESS FOR MAX SLOWNESS AND MIN SIGL=* 
1,F15.5) 

FORMAT (1H ,"MAX THICKNESS OF LVL!,15, *=',F15.5) 

FORMAT (1H ,*MAXIMUM THICKNESS FOR MAX SLOWNESS IS? 
1,F15.5) 

FORMAT(1H ,*MAX SLOWNESS FOR MAX THICKNESS=",F15.10) 

FORMAT (1015) 

FORMAT (2F15. 3) 

FORMAT (3F20.0) 

FORMAT (1H ,2F10.3,/) 

FORMAT (2F 10.3) 

FORMAT(1H ,*NU,NL, PU, PL, AU, BU,CU,AL, BL,CL,SAM,NDIV 
1, VSUR, ILU, ILL, NREGS, NVELS, NLVZS, SIGU, SIGL, QU, QL 
1, UMAX, MUST*) 

FORMAT (1H ,3D20. 10) 

FORMAT(1H ,5D20.10,15) 

FORMAT (5F15.5,15) 

R=6371.0D0 

READ (5,100) NU,NL 

READ(5,101) (PU(I) ,I=1,NU) 

READ(5,101) (PL(I) ,1=1,NL) 

MU=NU-1 

ML=NL-1 

READ(5,102) ((AU(I) ,BU(I) ,CU(I)) ,I=1,MU) 

READ (5,102) ((AL(ZI),BL(I) CL (I)),I=1, ML) 

READ(5,100) (ILU(I) ,I=1, 40) 

READ (5,100) (ILL (I) ,I=1,ML) 

DO 2228 JV=1,NU 

PU (JV) =PL (JV) 

AU (JV) =AL (JV+1) 

CONTINUE 

DO 2229 JV=1,40 
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ILU (JV) =1 

ILL (JV) =1 

BU (JV) =0.0 

BL (JV) =0.0 

CU (JV) =0.0 

CL(JV) =0.0 
2229 CONTINUE 
READ (5,101) SAM 
READ(5,100) NDIV 
READ(5,101) VSUR 
READ(5,100) NREGS 
READ(5,100) (NVELS (I) ,I=1,NREGS) 
READ(5,100) NLVZS 
IF (NLVZS.EQ.0) GO TO 448 
READ(5,7000) ((SIGU (I) ,SIGL (I) ,QU (I) ,QL(1I) , UMAX (I) 
, MUST (I)) ,I=1,NLVZS) 
448 CONTINUE 

PD=PL (NL) 

WRITE (6, 4343) 

WRITE (6,100) NU,NL 

WRITE (6,4344) (PU(I) ,I=1,N0) 

WRITE (6,4344) (PL(I) ,I=1,NL) 
Y344 FORMAT(1H ,5D20.9) 
WRITE (6,4345) (AU(I) ,-BU(I) ,CU(I) ,1=1, MU) 
WRITE (6,4345) (AL(I),BL(I) ,CL(1I) ,I=1,ML) 
WRITE (6,4345) SAM | 
WRITE(6,100) NDIV 
WRITE (6,4345) VSUR 
WRITE(6,100) (ILU(I) ,I=1,MU) 
WRITE(6,100) (ILL(I),I=1,/ML) 
WRITE (6,100) NREGS, NVELS,NLVZS 
WRITE (6,4346) (SIGU(I),SIGL(I) ,QU(I) ,QL(I) , UMAX (I) 
oMUST(I)) ,1=1,NLVZS) 
IF (NLVZS.EQ.0) QU (1) =PL (NL) 
DIV=(1.0D0-QU (1) ) /(NVELS (1) #1) 
MMM=NVELS (1) 
DO 10 J=1,MMM 
A (J) =1.0D0-DIV«Jd 
B (J) =A (J) 
10 CONTINUE 

NV=NVELS (1) 

IF(NLVZS.EQ.0) GO TO 40 

II=0 

IF(NREGS.EQ-2) GO TO 30 

IJK=NREGS -1 

DO 20 I=2,1IdK 

IF (II.£Q.1) GO TO 60 

DIV= (QU (I/2-QL (I/2) ) / (NVELS (I-1) 

MMN=NVELS (T) 

DO 70 J=1,MMN 

A (J+NV) =QU (I/2- (J-1) «DIV 
70 B (J+NV) =A (J#NV) 

A (1+NV) =QU (1/2) 

B (1#NV) =QU (1/2) 
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A (NVELS (1) #NV) =QL (1/2) 
B (NVELS (I) #NV) =QL (1/2) 
GO TO 19 
60 DIV=(QL( (1-1) /2-QU (( (I-1) /2) #1) ) / (NVELS (I) #1) 
MNN=NVELS (I) 
DO 80 J=1,MNN 
A (J+NV) =QL ( (I-1) /2-d«DIV 
80 B(J#tNV)=A(J+4NV) 
19 NV=NV+NVELS (I) 
II=II-1 
IF (II.EQ.-1) II=1 
20 CONTINUE 
30 CONTINUE 
IF(II.£Q.0) GO TO 90 
DIV= (QL ( (NREGS-1) /2-PD) / (NVELS (NREGS) +1) 
IXXX=NVELS (NREGS) 
DO 92 J=1,IXxX 
A (J+NV) =QL ( (NREGS-1) /2-J «DIV 
92 B(J+tNV)=A(J+NV) 
NV=NV+NVELS (NREGS) 
GO TO 40. 
90 CONTINUE 
DIV= (QU (NREGS/2-PD) / (NVELS (NREGS) ) 
IXX=NVELS (NREGS) 
DO 93 J=1,1XxX 
A (J+NV) =QU (NREGS/2- (J-1) *DIV 
93 B(J#NV)=A(J+#NV) 
NV=NV4+NVELS (NREGS) 
40 CONTINUE | 
CALL LOWLIM(A(J) ,VA (J) ,BBEST (J) ) 
CALL UPPLIM(B(J) ,VB(J) ,ABEST (J) ) 
2000 CONTINUE 
IF(NREGS.EFQ.1) GO TO 4000 
MV1=NVELS (1) +1 
MV2=NVELS (1) +NVELS (2) 
NV1=2 
NV2=NVELS (1) #+NVELS (2) 
DO 800 K=1,NLVZS 
DO 79 J=NV1,NV2 
79 ZIF(VA(J)-LT.VA(J-1)) VA(J) =VA(J-1) 
JJ=0 
DO 78 J=MV1,MV2 
JIJ=JI+1 
78 SAVE (K,JJ) =VA (J) 
DO 77 J=MV1,MV2 
77  VA(J)=VA (MV1) 
IF(II.E£O0.0.AND.K.EQ.NLVZS) GO TO 800 
IF (II. £Q.1.AND.K.EQ.NLVZS) GO TO 3000 
LV1=NV2+1 
LV2=NV2+NVELS (2%K+1) +NVELS (2%K+ 2) 
DO 76 J=LV1,LV2 
IF (BBEST (J) .GT.QL(K)) GO TO 75 
SKE= (A (J) +BBEST (J) ) /2.0D0 
DO 55 KV=1,K 
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ARG1=DSQRT (OMAX (KV) *xUMAX (KV-OL (KV) «OL (KV) ) 
ARG2=DSQRT (QL (KV) *OL (KV-SKE*SKE) 
VA (J) =VA (J— (SIGU (KV) * (1. OD0- (ARG2/ARG1) *DATAN2 (ARG1 
1, ARG2))) *2.0D0/ (3.14159 26535897 9D0xARG2) 
55 CONTINUE 
GO TO 76 
75 VA(J)=VA(J-1) 
76 CONTINUE 
GO TO 799 
3000 CONTINUE 
LV1=NV2+1 
LV2=NV2+NVELS (2*K+1) 
DO 74 J=LV1,LV2 
IF (BBEST (J) .GT.QL(K)) GO TO 73 
SKE=(A(J) #BBEST (J) ) /2-0D0 
DO 99 KV=1,K 
ARG 1=DSQRT (UMAX (KV) «UMAX (KV-OL (KV) «QL (KV) ) 
ARG2=DSQRT (OL (KV) «OL (KV-SKE*SKE) - 
VA (J) =VA (J— (SIGU (KV) « (1. ODO— (ARG2/ARG1) xDATAN2 (ARG1,A 
1RG2) )) #2. 0D0/ (3. 14159 265358979D0«xARG2) 
99 CONTINUE 
GO TO 74 
73 +VA(d3)=VA(J-1) 
74 CONTINUE 
| DO 34 J=LV1,LV2 
34 IF(VA(J).LT.VA(J-1)) VA(J) =VA(J-1) 
GO TO 800 
799 CONTINUE 
IF (K.-EQ.NLVZS) GO TO 800 
MV1=MV2+NVELS (2«K+1) +1 
MV2=MV24#NVELS (2%K+1) #+NVELS (2%K+2) 
NV1I=NV2+1 
NV2=NV2¢NVELS (2*K+1) #*NVELS (24K+ 2) 
800 CONTINUE | 
MV1=NVELS (1) +1 
MV2=NVELS (1) #NVELS (2) 
NV1=MV2+1 
DO 900 K=1,NLVZS 
HMAX (K) =0.0D0 
JJ=0 
DO 5005 J=MV1,MV2 
JJ=JI+1 
TH=VB (J-SAVE(K,JJ) 
IF (TH. GT.HMAX (K)) HMAX(K) =TH 
5005 CONTINUE 
WRITE (6,41) K,HMAX (K) 
THM=SIGU (K) /DSORT (UMAX (K) *UMAX (K-QU (K) «QU (K) ) 
SLOW=DSQRT ( (SIGU (K) xSIGU (K) ) / (HMAX (K) #HMAX (K)) 
1+QU (K) *QU (K)) 
WRITE (6,42) THM 
WRITE (6,43) SLOW 
IF (MUST (K).£Q.0) GO TO 54 
THM=SIGL (K) /DSORT (UMAX (K) *UMAX (K-QU (K) *QU (K) ) 
SLOW=DSOQRT ( (SIGL (K) *SIGL (K)) / (HMAX (K) *HMAX (K)) 
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900 
4000 
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1+ QU (K) «QU (K)) 

1U (K)) 

WRITE(6,39) THM 

WRITE (6,38) SLOW 

UMIN (K) =DSQRT ( (SIGL (K) *«SIGL (K)) / (HMAX (K) xHMAX (K) ) 
1+Q0 (K) «QU (K) ) 

WRITE(6,37) UMIN (K) 
IF(IT.EQ.0.AND.K.-EQ.NLVZS) GO TO 900 
DO 88. J=NV1,NV 

SKE= (ABEST (J) +B (J) ) /2.0D0 
ARG1=DSQRT (OMIN (K) x UMIN (K-QD (K) «QU (K) ) 
ARG2=DSQRT (QU (K) «QU (K-SKE*SKE) 

VB(J) =VB(J- (SIGL (K) « (1.0D0—-(ARG2/ARG1) xDATAN2 (ARG1,AR 
1G2))) *2.0D0/ (3.14159 26535897 9x ARG2) 
CONTINUE 

IF (11. £EQ.1.AND.K.EQ.NLVZS) GO TO 900 
CONTINUE 
MV1=MV24NVELS (2% K4#1) +1 
MV2=MV2+NVELS (2%K+1) +NVELS (2%K+ 2) 
NV1I=MV2+1 

CONTINUE 

CONTINUE 


-DO 3333 J=2,NV 


IF(VA(J)..LT.VA(J-1)) VA(J) =VA(J-1) 
CONTINUE 
K=NV 

J2=NV-1 

DO 3334 J=1,J92 

K=K-1 

IF (VB(K) .GT.VB(K+1)) VB(K) =VB(K+1) 
DO 3335 J=1,NV 

FACT=DEXP (—VSUR«VA (J) /R) 
VEL=VSURxFACT/A(J) 

DEP=Rx (1.0D0-FACT) 

WRITE (6,2222) VEL,DEP 
WRITE (7,2223) VEL,DEP 
PACT=DEXP (—VSUR«VB (J) /R) 
VEL=VSUR*FACT/B (J) 

DEP=Rx (1.0D0-FACT) 

WRITE (6,2222) VEL,DEP 

WRITE (7,2223) VEL,DEP 
CONTINUE 

CALL EXIT 

END 
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SUBROUTINE LOWLIM(A, VBEST, BBEST) 
DOUBLE PRECISION A, VBEST, BBEST, SAMP,SAM,B1, B2 
DOUBLE PRECISION VAL1,VAL2,BOLD,OLD,SDIV,B, VAL 
DOUBLE PRECISION PU (210) ,PL (210) 
DOUBLE PRECISION AU (210) ,BU (210) ,CU (210) 
DOUBLE PRECISION AL(210),BL (210) ,CL (210) 
DIMENSION ILL (210) , ILU (210) 
COMMON PL, PU 
COMMON AU,BU ,CUy Aly BL, CL 
COMMON SAM 
COMMON NDIV,NL,NU 
COMMON ILL,ILU 

1000 FORMAT (1H ,3F15.5) 

DO 10 II=2,NL 
T1I7=I1I 
IF(A.GE.PL(II)) GO TO 11 

10 CONTINUE 

11. NRLA=I17-1 
SAMP=SAM 

1 B1=A4SAMP 
B2=A+2.0D0*SAMP 
IF (B2.LE.1.0D0) GO TO 2 
3 SAMP=SAMP/2.0D0 
GO TO 1 
2 CALL LLIM(A,NRLA,B1,VAL1) 
CALL LLIM(A, NRLA, B2, VAL2) 
IF(VAL2.GT.VAL1) GO TO 4 
GO TO 3 
4 BOLD=B1 
OLD=VAL1 
VAL1=VAL2 
B1=B2 
B2=B2+SAMP 
IF(B2.-LT.1.0D0) GO TO 5 
B2=1.0D0 
CALL LLIM(A, NRLA,B2, VAL2) 
IF(VAL2.LE.VAL1) GO TO 6 
BBEST=B2 
VBEST=VAL2 
GO TO 7 
5S CALL LLIM(A,NRLA, B2,VAL2) 
IF(VAL2.GT.VAL1) GO TO 4 
6  VBEST=VAL1 
BBEST=B1 
SDIV=(B2-BOLD) / (NDIV+1) 
DO 90 K=1,NDIV 
B=BOLD+KxSDIV 
CALL LLIM(A,NRLA,B, VAL) 
IF(VAL.GT.VBEST) BBEST=B 
IF(VAL.GT.VBEST) VBEST=VAL 
90 CONTINUE 
7 CONTINUE 
WRITE (6, 1000) A,VBEST,BBEST 
RETURN 
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SUBROUTINE LLIM(A,NRLA,B, V) 

DOUBLE PRECISION A,B,V,P1,P2,P3,P4,R1,R2, R3, RY 
DOUBLE PRECISION 01,02,03,04,T3,T4,239Z4,V1,V2 
DOUBLE PRECISION SAM 

DOUBLE PRECISION PU (210) ,PL (210) 

DOUBLE PRECISION AU (210), BU (210) ,CU (210) 
DOUBLE PRECISION AL (210) , BL (210) ,CL (210) 
DIMENSION ILL (210) , ILU (210) 

COMMON PL, PU 

COMMON AU,BU,CU,AL,BL,CL 

COMMON SAM 

COMMON NDIV,NL,NU 

COMMON ILL, ILU 

DO 10 I=2,NL 

I7=I 

IF(B.GT.PL(I)) GO TO 11 

CONTINUE 

NRBF=1I7-1 

DO 20 I=2,NU 

I7=1 

IF(B.GE.PU(I)) GO TO 21 

CONTINUE 

NRBS=I7-1 

IR=NRLA-NRBF 

P2=B 

IF(IR.GT.0) P2=PL(NRLA) 

R2=A/P2 

Q2=P2*P2-AxA 

V1=AL(NRLA) *DARCOS (R2) #AxBL (NRLA) * (DLOG(P2+DSORT (02) ) 


1-DLOG{A)) 


IF (ILL (NRLA) -EQ.1) GO TO 110 
V1I=V1I+A«xCL (NRLA) «DSQRT (Q2) 
CONTINUE 

IF(IR.EQ-0) GO TO 40 

P1=PL (NRBF+1) 

P2=B 

R2=A/P2 

Q2=P2*P2-AxA 

Q1=P1*xP1-AxA 

R1=A/P1 

V1=V14AL (NRBF) * (DARCOS (R2-DARCOS (R1) ) #AxBL (NRBF) *DLOG 


1( (P2+DSQRT (Q2)) /(P1*+DSQRT (Q1))) 


IF (ILL (NRBF).EQ.1) GO TO 111 


~V1=V14A«xCL (NRBF) « (DSQRT (Q2-DSQRT (Q1) ) 


CONTINUE 

IF (IR.EQ.1) GO TO 40 
NT=IR-1 

DO 50 J=1,NT 

P1=PL (NRLA-J+1) 
P2=PL (NRLA-J) 
R1=A/P 1 

R2=A/P2 

Q1=P1*P1-AxA 
Q2=P2*P2-AxA 
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V1=V1+AL (NRLA-J) * (DARCOS (R2-—DARCOS (R1)) +AxBL (NRLA-J) x 
IDLOG ( (P2+DSQRT (Q2)) / (P1#DSQRT(Q1))) 

IF (ILL (NRLA-J) .EQ.1) GO TO 112 

V1=V 1+AxCL (NRLA-J) « (DSQRT (Q2—DSQRT (Q1) ) 

CONTINUE 

CONTINUE 

CONTINUE 

V¥2=0.0D0 

IF (B.EQ.1.0D0) GO TO 60 

P3=B 

P4=PU (NRBS) 

R4=A/P4 

R3=A/P3 

Q3=P3*P3—-AxA 

Q4U=P4xPU—-AxA 

TU=B/P4 

Z4=P4UxP4U-BxB 

V2=AU (NRBS) « (DARCOS (R4U—DARCOS (R3-—DARCOS (TY) ) +BU (NRBS 
1) « (Ax (DLOG (P4#DSORT (Q4) —-DLOG (P3*+DSQRT (Q3) ) ) + Bx (DLOG (P 
13-DLOG (P4+DSQORT (Z4) ))) 

IF (ILU (NRBS) -EQ.1) GO TO 113 

V2=V2+CU (NRBS) «(Ax (DSQORT (Q4-DSORT (03) -BxDSOQRT (Z4) ) 
CONTINUE 

IF (NRBS.EQ.1) GO TO 60 

MT=NRBS-1 

DO 70 J=1,MT 

P3=PU (NRBS+1-J) 

P4=PU (NRBS-J) 

R4=A/P4 

R3=A/P3 

T3=B/P3 

T4U=B/P4 

Q4U=P4UxP4-AxA 

Q3=P3xP3-AxA 

Z3=P3«P3-Bx«xB 

Z4=P4xP4-BxB 

V2=V2+AU (NRBS-J) x (DARCOS (R4—DARCOS (R3-DARCOS (T4) +DAR 
1COS (T3) ) +BU (NRBS-J) x (Ax (DLOG (P4*#DSORT (Q4) -DLOG (P3+DSQ 
1RT(03))) +Be (DLOG (P3+DSQRT (23) -DLOG (P44+DSQRT (24) ))) 
IF (ILU (NRBS-J) .EQ.1) GO TO 114 

V2=V2+CU (NRBS-J) « (Ax (DSQRT (Q4U-—DSORT (03) —Bx (DSORT (Z4) 
1-DSORT (Z3))) 

CONTINUE 

CONTINUE 

CONTINUE 

V=(V1+V2) «2. 0D0/ (B-A) 

V=V/3. 1415926535897 9D0 

RETURN 

END 
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SUBROUTINE UPPLIM(B,VBEST,ABEST) 


DOUBLE 
DOUBLE 
DOUBLE 
DOUBLE 
DOUBLE 


PRECISION B,VBEST, ABEST, SAMP,SAM,A1,A2,PD 
PRECISION VAL1,VAL2,AOLD,OLD,SDIV,A, VAL 
PRECISION PU (210) , PL (210) 

PRECISION AU(210) , BU (210) ,CU (210) 
PRECISION AL (210) ,BL(210) ,CL (210) 


DIMENSION ILL (210) , ILU (210) 


COMMON 
COMMON 
COMMON 
COMMON 
COMMON 


PL, PU 

AU, BU,CU, AL, BL,CL 
SAM 

NDIV,NL,NU 

ILL, ILU 


FORMAT (1H ,3F15.5) 
PD=PU (NU) 
DO 10 II=2,NU 


II7=1I1 


IF(B.GT.PU(II)) GO TO 11 
CONTINUE 

NRUBF=11I7-1 

DO 20 JJ=2, NL 


JIIT=IJI 


IF(B.GE.PL(JJ)) GO TO 21 
CONTINUE 

NRUBS=JJ7-1 

SAMP=SAM 

A1=B-SAMP 


A2=B-2 @ 


ODO*xSAMP 


IF(A2.GE.PD) GO TO 2 
S AMP=SAMP/2.0D0 


GO TO 1 


CALL ULIM(B, NRUBF,NRUBS,A1,VAL1) 
CALL ULIM(B, NRUBF, NRUBS, A2,VAL2) 
IF(VAL1.GT.VAL2) GO TO 4 


GO TO 3 
AOLD=A1 


OLD=VAL1 
VAL1=VAL2 


Al=A2 


A2=A2-SAMP 
IF(A2.-GT.PD) GO TO 5 


A2=PD 


CALL ULIM(B,NRUBF,NRUBS,A2,VAL2) 
IF(VAL2.GE.VAL1) GO TO 6 
ABEST=PD 

VBEST=VAL2 


GO TO 7 


CALL ULIM(B, NRUBF, NRUBS, A2, VAL2) 
IF(VAL2.LT.VAL1) GO TO 4 
VBEST=VAL1 

ABEST=A1 

SDIV= (AOLD-A2) / (NDIV+#1) 

DO 90 K=1,NDIV 


A=AOLD- 


K*SDIV 


CALL ULIM(B,NRUBF, NRUBS, A, VAL) 
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IF(VAL.LT.VBEST) ABEST=A 
IF(VAL.LT.VBEST) VBEST=VAL 
CONTINUE 

CONTINUE 

WRITE(6,1000) B,VBEST,ABEST 
RETURN 

END 
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SUBROUTINE ULIM(B,NRUBF, NRUBS,A,V) 
DOUBLE PRECISION B,A,V,P1,P2,P3,P4,R1,R2,R3, RU 
DOUBLE PRECISION Q1,02,03,04,T3,T4,23,24 
DOUBLE PRECISION V1,V2 

DOUBLE PRECISION SAM 

DOUBLE PRECISION PU (210) , PL (210) 

DOUBLE PRECISION AU (210) , BU (210) ,CU (210) 
DOUBLE PRECISION AL (210) ,BL(210) ,CL (210) 
DIMENSION ILL(210),ILU (210) 

COMMON PL,PU 

COMMON AU, BU,CU, AL, BL,CL 

COMMON SAM 

COMMON NDIV,NL,NU 

COMMON ILL,ILV 

DO 10 I=2,NU 

I7=1 

IF(A.GE.P0 (I)) GO TO 11 

CONTINUE 

NRA=I7-1 

TR=NRA-NRUBF 

P2=B 

IF(IR.GT.0) P2=PU(NRA) 

R2=A/P2 

Q2=P2*P2-AxA 

V1=AU (NRA) xDARCOS (R2) +AxBU (NRA) * (DLOG (P2*DSQRT (Q2) -DL 
10G (A) ) 

IF (ILU (NRA) -EQ.1) GO TO 110 
V1=V1+AxCU (NRA) «DSQRT (Q2) 

CONTINUE 

IF(IR.~EQ.0) GO TO 40 

P1=PU (NRUBF+t1) 

P2=B 

R1=A/P1 

R2=A/P2 

Q2=P2*P2—-AxA 

Q1=P1*xP1-AxA 

V1=V1+AU (NRUBF) x (DARCOS (R2-DARCOS (R1) ) #AxBU (NRUBF) «DL 
10G ((P2+DSQRT (Q2)) / (P1*+DSQRT (Q1) ) ) 

IF(ILU (NRUBF).EQ.1) GO TO 111 
V1=V1+AxCU (NRUBF) « (DSQRT (Q2-DSQRT (Q1) ) 
CONTINUE 

IF(IR.EQ.1) GO TO 40 

NT=IR-1 

DO 50 J=1,NT 

P1=PU (NRA-J+1) 

P2=PU (NRA-J) 

R1=A/P1 

R2=A/P2 

Q1=P1*xP1—-AxA 

Q2=P2*xP2—-AxA 

V 1=V1+AU (NRA-J) x (DARCOS (R2-DARCOS (R1) ) tAxBU (NRA-J) «DL 
10G((P2+DSORT (02) ) / (P1*DSQRT (Q1) )) 

IF (ILU (NRA-J) -EQ.1) GO TO 112 
V1=V1+A«xCU (NRA-J) «x (DSQRT (Q2-DSQRT (Q1) ) 
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CONTINUE 

CONTINUE 

CONTINUE 

V2=0.0D0 

IF (B.EQ.-1.0D0) GO TO 60 

P3=B 

P4=PL (NRUBS) 

R3=A/P3 

R4=A/P4 

T4=B/P4 

Q3=P3*%P3-AxA 

Q4U=P4xP4-AxA 

Z4=P4x%P4U-BxB 

V2=AL (NRUBS) « (DARCOS (R4-—DARCOS (R3—DARCOS (T4) ) +BL (NRU 
1BS) « (Ax (DLOG (P4+DSQRT (04) -DLOG (P3+DSQRT (Q3) ) ) + Bx (DLOG 
1 (P3-DLOG (P4+DSQRT (Z4) ))) 

IF (ILL (NRUBS) .EQ.1) GO TO 113 

V2=V2+CL (NRUBS) x (Ax (DSORT (QU-DSQRT (03) —-BxDSQRT (Z4) ) 
CONTINUE 

IF (NRUBS.EQ.1) GO TO 60 

MT=N RUBS-1 

DO 70 J=1,MT 

P3=PL (NRUBS+ 1-J) 

P4=PL (NRUBS-J) 

R4=A/P4 

R3=A/P3 

T3=B/P3 

T4=B/P4 

Q3=P3*xP3-AxA 

Q4U=P4xP4—-AxA 

Z3=P3%P3-BxB 

Z4=P4xP4—-BxB 

V2=V2+AL (NRUBS-J) « (DARCOS (R4—DARCOS (R3-DARCOS (T4) +DA 
1RCOS (T3) ) #BL (NRUBS—J) « (Ax (DLOG (P4+#DSQRT (Q4) —-DLOG (P3+ 
IDSORT (03) )) +Bx (DLOG (P3+DSQRT (Z3) -DLOG (P4#DSQRT (Z4)))) 
IF (ILL (NRUBS-J) .£Q.1) GO TO 114 

V2=V2+CL (NRUBS-J) «(Ax (DSORT (Q4—-DSORT (03) —Bx (DSORT (24 
1-DSQRT (Z3))) 

CONTINUE 

CONTINUE 

CONTINUE 

V=(V1+V2) *2.0D0/ (B-A) 

V=V/3. 14159265358979D0 

RETURN 

END 
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